Review of Monte Carlo modeling of light transport in tissues

Abstract. A general survey is provided on the capability of Monte Carlo (MC) modeling in tissue optics while paying special attention to the recent progress in the development of methods for speeding up MC simulations. The principles of MC modeling for the simulation of light transport in tissues, which includes the general procedure of tracking an individual photon packet, common light–tissue interactions that can be simulated, frequently used tissue models, common contact/noncontact illumination and detection setups, and the treatment of time-resolved and frequency-domain optical measurements, are briefly described to help interested readers achieve a quick start. Following that, a variety of methods for speeding up MC simulations, which includes scaling methods, perturbation methods, hybrid methods, variance reduction techniques, parallel computation, and special methods for fluorescence simulations, as well as their respective advantages and disadvantages are discussed. Then the applications of MC methods in tissue optics, laser Doppler flowmetry, photodynamic therapy, optical coherence tomography, and diffuse optical tomography are briefly surveyed. Finally, the potential directions for the future development of the MC method in tissue optics are discussed.


Introduction
Monte Carlo (MC) methods are a category of computational methods that involve the random sampling of a physical quantity. 1,2 The term "the Monte Carlo method" can be traced back to 1940s, 1 in which it was proposed to investigate neutron transport through various materials. Such a problem cannot be solved by conventional and deterministic mathematical methods. Due to its versatility, this method has found applications in many different fields 3 including tissue optics. It has become a popular tool for simulating light transport in tissues for more than two decades 4 because it provides a flexible and rigorous solution to the problem of light transport in turbid media with complex structure. The MC method is able to solve radiative transport equation (RTE) with any desired accuracy, 5 assuming that the required computational load is affordable. For this reason, this method is viewed as the gold standard method to model light transport in tissues, results from which are frequently used as reference to validate other less rigorous methods such as diffuse approximation to the RTE. 6,7 Due to its flexibility and recent advances in speed, the MC method has been explored in tissue optics to solve both the forward and inverse problems. In the forward problem, light distribution is simulated for given optical properties, whereas in the inverse problem, optical properties are estimated by fitting the light distribution simulated by the MC method to experimentally measured values.
In this review paper, the principles of MC modeling for the simulation of light transport in tissues, including the general procedure of tracking an individual photon packet, common light-tissue interactions that can be simulated such as light absorption and scattering, frequently used tissue models, common contact and noncontact illumination and detection setups, and the treatment of time-resolved and frequency-domain optical measurements, are described in detail to help interested readers achieve a quick start. Following that, a variety of methods for speeding up MC simulations, including scaling methods, perturbation methods, hybrid methods, variation reduction techniques, parallel computation, and special methods for fluorescence simulations, and their respective advantages and disadvantages are discussed. Then the biomedical applications of MC methods, including the simulation of optical spectra, estimation of optical properties, simulation of optical measurements in laser Doppler flowmetry (LDF), simulation of light dosage in photodynamic therapy (PDT), simulation of signal source in optical coherence tomography (OCT) and diffuse optical tomography (DOT), are surveyed. Finally, the potential directions for the future development of MC methods are discussed, which are based on their current status in the literature survey and the authors' anticipation. It should be pointed out that this review is intended to give a general survey on the capability of MC modeling in tissue optics while paying special attention on methods for speeding up MC simulations since the timeconsuming nature of common MC simulations could limit its applications. In the general procedure of MC modeling, light transport in tissues is simulated by tracing the random walk steps that each photon packet takes when it travels inside a tissue model. For each launched photon packet, an initial weight is assigned as it enters the tissue model, as illustrated in Fig. 1. The step size will be sampled randomly based on the optical properties of the tissue model. If it is about to hit a boundary, either of the following two methods could be used to handle this situation. In the first method, the photon packet will either transmit through or be reflected from the boundary. In the second method, a fraction of the photon packet's weight will always be reflected and the remaining fraction of the photon packet's weight will transmit through. The probabilities of transmission or reflection in the first method, and the fraction of the photon packet's weight transmitting through or being reflected in the second method, are governed by Snell's law and Fresnel's equations. At the end of each step, the photon packet's weight is reduced according to the absorption probability; meanwhile, the new step size and scattering angle for the next step will be sampled randomly based on their respective probability distributions. The photon packet propagates in the tissue model step by step until it exits the tissue model or is completely absorbed. Once a sufficient number of photon packets are launched, the cumulative distribution of all photon paths would provide an accurate approximation to the true solution of the light transport problem and the contribution averaged from all photons can be used to estimate the physical quantities of interest.

Common Light-Tissue Interactions in MC Modeling
Several types of common light-tissue interactions, including light absorption, elastic scattering, fluorescence and Raman scattering, have been simulated by the MC methods previously. The absorption coefficient μ a (unit: cm −1 ) and the scattering coefficient μ s (unit: cm −1 ) are used to describe the probability of absorption and scattering, respectively, occurring in a unit path length. The anisotropy factor g, which is defined as the average cosine of scattering angles, determines the probability distribution of scattering angles to the first-order approximation. In addition, the refractive index mismatch between any two regions in the tissue model or at the air-tissue interface will determine the angle of refraction. The fraction of photon packet weight that, after traveling in the medium, escapes from the same side of the tissue model as the incident light is scored as diffuse reflectance. In contrast, the fraction of photon packet weight that travels through the medium and escapes from the other side of the tissue model is scored as transmittance. 5,8,9 To simulate fluorescence emission, one additional parameter, which is fluorescence quantum yield, 10,11 needs to be incorporated to describe the probability that the absorbed photon packet weight can be converted to a fluorescence photon at a different wavelength. If time-resolved fluorescence is simulated, the lifetime of fluorescence needs to be defined. The initial direction of the fluorescence photon is isotropic due to the nature of fluorescence emission. As illustrated in Fig. 2, the MC modeling of fluorescence propagation in tissues involves three steps. [11][12][13] The first step involves a general MC simulation to simulate light propagation with optical properties at the excitation wavelength. In the second step, a fluorescence photon may then be generated upon the absorption of an excitation photon with a probability defined by the quantum yield and time delay defined by the lifetime of fluorescence. The third step again involves a general MC simulation to simulate fluorescence light propagation with optical properties at the emission wavelength. It is clear that simulated fluorescence from a tissue model will be related to the absorption and scattering properties of the tissue model in addition to the fluorescence quantum yield and lifetime. Fluorescence simulation is typically much more time-consuming than the simulation of diffuse reflectance due to extra fluorescence photon propagation.
To simulate Raman emission, a parameter similar to fluorescence quantum yield, named as Raman cross-section, [14][15][16][17] is needed to describe the probability that a Raman photon will be generated at each step. A phase function for Raman photons needs to be determined. The MC simulation procedure for Raman light propagation will be similar to that for fluorescence.
Bioluminescence refers to the phenomenon of living creatures producing light, which results from the conversion of chemical energy to bioluminescence photons 18 and which has  also been investigated using MC modeling. Because bioluminescence does not need an external light source for excitation, the first step in MC simulation of bioluminescence is to generate bioluminescence photons package according to the distribution of bioluminescence sources. 19,20 After that, the simulation of bioluminescence photon propagation in a tissue model is exactly the same as the simulation of diffuse reflectance.
If the polarization property of light is considered in MC modeling, the polarization of a photon can be represented by Stokes vectors and the polarimetry properties of the tissue model can be described by Jones matrix or Mueller matrices, 21,22 which will not be expanded in this review.

Common Tissue Models in MC Modeling
Common tissue models used in MC simulations include the homogeneous and nonhomogeneous tissue models. The optical properties in a homogeneous tissue model are equal everywhere. 4,23 In contrast, the optical properties in a nonhomogeneous tissue model vary with the tissue region. The following survey is focused on nonhomogeneous tissue models because of its high preclinical and clinical relevance.
The most commonly used nonhomogeneous tissue model is perhaps the multilayered tissue model, 8,9,12,13,[24][25][26][27][28][29][30] which is frequently employed to mimic epithelial tissues. In a multilayered tissue model, each tissue layer is assumed to be flat with uniform optical properties and it is infinitely large on the lateral dimension. This assumption works fine when the source-detector separation is small so that the spatial variation in the optical properties within the separation is negligible. However, it could cause significant errors if the optical properties change significantly in a small area, such as in dysplasia or early cancer 31 and port wine stain (PWS) model. 32 To overcome this limitation, tissue models including heterogeneities with well-defined shapes have been used to mimic complex tissue structures from different organs. For example, Smithies et al. 32 and Lucassen et al. 33 independently proposed MC models in which simple geometric shapes were incorporated into layered structures to model light transport in PWS model. In their PWS models, infinitely long cylinders were buried in the bottom dermal layer to mimic blood vessels. Wang et al. 34 reported an MC model in which a sphere was buried inside a slab to model light transport in human tumors. Zhu et al. 31,35 proposed an MC model in which cuboid tumors were incorporated into layered tissues to model light transport in early epithelial cancer models including both squamous cell carcinoma and basal cell carcinoma.
Voxelated tissue models have been also explored to simulate irregular structures. Pfefer et al. 36 reported a three-dimensional (3-D) MC model based on modular adaptable grids to model light propagation in geometrically complex biological tissues and validated the code in a PWS model. Boas et al. 37 proposed a voxel-based 3-D MC model to model arbitrary complex tissue structures and tested the code in an adult head model. Patwardhan et al. 38 also proposed a voxel-based 3-D MC code for simulating light transport in nonhomogeneous tissue structures and tested the code in a skin lesion model. The three voxelbased MC codes above showed great flexibility in a range of applications. However, to model tissue media with curved boundaries in a voxel-based MC model, the grid density will have to be increased, which requires more memory and computation. A few other approaches have been explored to accommodate this situation. Li et al. 19 reported a public MC domain, named mouse optical simulation environment, to model bioluminescent light transport in a living mouse model. The mouse model consists of several segmented regions that are extended from several building blocks such as ellipses, cylinders, and polyhedrons. This platform is particularly suitable for small animal imaging. Margallo-Balbas et al. 39 and Ren et al. 40 have developed triangular-surface-based MC methods to model light transport in complex tissue structures. The triangular-surface-based approach allows an improved approximation to the interfaces between domains, but it is not able to model complex media with continuously varying optical properties. Moreover, it could be time-consuming to determine ray-surface intersection because a range of triangles will have to be scanned. To overcome the limitations associated with triangular-surface-based MC method, most recently, Shen et al. 41 as well as Fang 42 have presented mesh-based MC methods, by which one can model much more complex structures and situations.

Common Illumination and Detection Setups in MC Modeling
One important advantage of MC modeling, as compared to other non-numerical methods such as diffuse approximation, is its capability to faithfully simulate a variety of contact and noncontact illumination and detection setups for optical measurements. Note that the contact setup requires the direct contact between the tip of an optical probe and tissue samples. In contrast, the noncontact setup enables optical measurements from a tissue sample without directly contacting it.

Contact setup for illumination and detection
Fiber-optic probes are commonly used in contact illumination and detection configurations as demonstrated in many previous reports. 43 In general, these fiber-optic probes could be divided into two groups. In the first group, the same fiber or fiber bundle is used for both illumination and detection, 30,44,45 while in the second group, separate fibers are used for illumination and detection. 12,[46][47][48] There is no difference in the treatment of these two groups of fiber-optic probes from the point of view of modeling because the first group of probes can be viewed as two separate and identical fibers or fiber bundles for illumination and detection that happen to locate at the same spatial position.
The key parameters in simulated fiber-optic probes include the radii, numerical apertures (NA), tilt angles of illumination and detection fibers, and the center-to-center distance between the two sets of fibers (which is called the source-detector separation), as well as the refractive indices of these fibers relative to that of the tissue model. The radius and NA of the illumination fiber in combination with the radial and angular distributions of photons coming out of the fiber define the locations and the incident angles of incident photons. For a commonly used multimode fiber, the spatial locations and incident angles of launched photons are typically assumed to follow uniform distribution and Gaussian distribution. Both spatial locations and incident angles need to undergo spatial coordinate transformation when the tilt angle of the illumination fiber is larger than zero. Here the tilt angle of a fiber refers to the angle of the fiber axis relative to the normal axis of the tissue model. The incident beam could be also assumed to be collimated or focused.
Light detection by a fiber usually contains two steps. The first step is to determine whether an exiting photon could enter the area defined by the radius of the detection fiber. If it is true, the second step is to determine whether the exiting direction of the photon falls within the acceptable angle of the detection fiber calculated from the NA and refractive index of the fiber. If the tilt angle of the detection fiber is larger than zero, the exiting location and angle are subject to spatial coordinate transformation.

Noncontact setup for illumination and detection
Noncontact setups usually employ various lenses for illumination and detection. In these setups, an adjunct lens or a combination of lenses is usually placed between a fiber-optic probe and the tissue sample to achieve noncontact measurements while maintaining the well-defined illumination and detection geometry. Jaillon et al. 49 proposed a method to simulate a beveled fiber-optic probe coupled with a ball lens to achieve depthsensitive fluorescence measurements from layered tissue models. Later, the same group incorporated a half-ball lens into the beveled fiber-optic probe to achieve the same purpose with a higher sensitivity. 50 Zhu et al. 51 proposed a method to simulate a fiber-optic probe coupled with convex lenses to achieve noncontact depth-sensitive diffuse reflectance measurements from early tumors in an epithelial tissue model. By manipulating the lens combination, an ordinary cone configuration and a special cone shell configuration were investigated. It was found that the cone shell configuration provides higher depth sensitivity to the tumor than the cone configuration.

Time-Resolved and Frequency-Domain MC Modeling
Time-resolved optical measurements such as fluorescence lifetime imaging (FLIM) 52 37,[57][58][59][60] The refractive index in each tissue region will influence the time that photons take to travel through. In the simulation of FLIM, it needs to be pointed out that the time delay from photon absorption to fluorescence generation should follow the probability density distribution defined by the fluorescence lifetime. 66,67 In the frequency-domain measurements, the modulation and/or phase delay of detected waves were analyzed. The modulation and phase delay can be simulated in either a direct approach 64 or an indirect approach, i.e., using Fourier transformation from a time-domain MC simulation. 65

Methods for the Acceleration of MC Simulation
While the MC method is the gold standard method to model light transport in turbid media, the major drawback of the MC method is the requirement of intensive computation to achieve results with desirable accuracy due to the stochastic nature of MC simulations, which makes it extremely time-consuming compared to other analytical or empirical methods. Significant efforts have been made to speed up the MC simulation of light transport in tissues during the past decades. These acceleration methods can be roughly divided into several categories as follows.

Scaling Methods
A typical scaling method requires a single or a few baseline MC simulations, in which the histories of survival photons such as trajectories or step sizes are recorded. Then, diffuse reflectance or transmittance for a tissue model with different optical properties can be estimated by applying scaling relations on the recorded photon histories. These methods take advantage of the fact that the scattering properties determine photon paths and the absorption property only influences the weights of survival photons. Graaff et al. 68 proposed a limited scalable MC method for fast calculation of total reflectance and transmittance from slab geometries with different optical properties. It was demonstrated that the trajectory information obtained in a reference MC simulation with a known albedo, i.e., μ s ∕ðμ a þ μ s Þ, can be used to find the total reflectance and total transmittance from slabs with other albedos. Kienle et al. 69 extended Graaff's theory to simulate space-and time-resolved diffuse reflectance from a semi-infinite homogeneous tissue model with arbitrary optical properties. Their approach was based on scaling (for different scattering coefficients) and re-weighting (for different absorption coefficients) a discrete representation of the diffuse reflectance from one baseline MC simulation in a nonabsorbing semi-infinite medium. It is powerful, but both the discrete representation and interpolation could introduce errors that are often amplified in scaling. Pifferi et al. 70 proposed a similar approach to estimate space-and time-resolved diffuse reflectance and transmittance from a semi-infinite homogeneous tissue model with arbitrary optical properties. Different from Kienle's method, the evaluation of reflectance and transmittance in Pifferi's approach is based on interpolation of results from MC simulations for a range of different scattering coefficients, and scaling is performed for absorption coefficients. This approach increases the accuracy of results for different scattering coefficients at the cost of a significantly increased number of baseline MC simulations. The methods reviewed above are fast, but the binning and interpolation involved introduce errors. In order to improve the accuracy of these methods, Alerstam et al. 60 improved Kienle's method by applying scaling to individual photons. In this method, the radial position of the exiting location and the total path length of each detected photon are recorded and the trajectory information of each photon will be individually processed to find the survival photon weight for tissue media with other sets of optical properties. Martinelli et al. 71 derived a few scaling relationships from the RTE, and their derivation showed that a rigorous application of the scaling method requires rescaling to be performed for each photon's biography individually. Two basic relations for scaling a survival photon's exit radial position r and exit weight w are listed in Eqs. (1) and (2) below. 47 in which r, w, μ t , and α are the exit radial position, exit weight, transport coefficients, and albedo in the baseline simulation, while r 0, w 0 , μ 0 t , and α 0 are those in the new simulations. N is the number of collisions recorded in the baseline simulation before the photon exits. Two relations essentially assume that the same set of random numbers sampled in the baseline simulation are also used in the new simulation and everything remains unchanged in two simulations, except the absorption and scattering coefficients.
Illumination and detection geometries have also been incorporated into the scaling procedure. Palmer et al. 47 extended Graaff's scaling method from illumination by a pencil beam to that by an optical fiber, and they also extended the original scaling method from the total reflectance to the reflectance detected by an optical fiber by combing scaling and convolution. Wang et al. 72 proposed two convolution formulas for the scaling MC method to calculate diffuse reflectance from a semi-infinite medium for a single illumination-detection fiber. Nearly all the previous papers about scaling dealt only with a homogeneous tissue model. Liu et al. 73 developed a method that applies the scaling method to multilayered tissue models. In this method, the homogeneous tissue model in a single baseline MC simulation is divided into multiple thin pseudo layers. The horizontal offset and the number of collisions that each survival photon experienced in each pseudo layer are recorded and used later to scale for the exit distance and exit weight of the photon in a multilayered tissue model with different set of optical properties. The method has been validated on both two-layered and three-layered epithelial tissue models.

Perturbation MC Methods
Similar to the scaling method, the perturbation MC (pMC) method requires one baseline simulation in which the optical properties are supposed to be close to the optical properties in the new tissue model so that the approximation made by perturbation is valid. 74 The trajectory information including the exit weight, path length, and number of collisions of each detected photon spent in the region of interest will be recorded in the baseline simulation. Then the relation between the survival weight in the baseline simulation and that in the new tissue model based on the perturbation theory, 75,76 i.e., is used to estimate diffuse reflectance from the tissue model in which the optical properties of the interesting region are perturbed. In Eq. (3), w, μ s , and μ t are the exit weight, scattering coefficient, and transport coefficient in the baseline simulation, while w new , μ 0 s , and μ 0 t are those in the new simulation. S and j are the photon path length and the number of collisions that a detected photon experienced in the perturbed region, respectively, recorded in the baseline simulation. It should be pointed out that the pMC is an approximation in nature, so its accuracy depends on the magnitude of difference in the optical properties between the perturbed optical properties in the new tissue model and the original optical properties in the baseline simulation.
In contrast, the scaling method is precise in nature regardless of the differences in optical properties because no approximation is made in scaling. One important advantage of the pMC is its simplicity and fast speed when the perturbed region is small, therefore it has been explored in the inverse problem of light transport to estimate optical properties in the perturbed region as surveyed below.
Sassaroli et al. 75 proposed two perturbation relations to estimate the temporal response in diffuse reflectance from a medium, in which scattering or absorbing inhomogeneities are introduced, from the trajectory information obtained from the baseline simulation of a homogeneous medium. Hayakawa et al. 76 demonstrated that the perturbation relation can be directly incorporated into a two-parameter Levenberg-Marquardt algorithm to solve the inverse photon migration problems in a two-layered tissue model rapidly. Recently, the same group 77 demonstrated the use of this method for extraction of optical properties in a layered phantom mimicking an epithelial tissue model for given experimental measurements of spatially resolved diffuse reflectance. This method was found effective over a broad range of absorption (50% to 400% relative to the baseline value) and scattering (70% to 130% relative to the baseline value) perturbations. However, this method requires both the thickness of the epithelial layer and the optical properties of one of the two layers.
Many other groups also proposed pMC-based methods for the recovery of the optical properties in various tissue models. Kumar et al. 78 have presented a pMC-based method for reconstructing the optical properties of a heterogonous tissue model with low scattering coefficients and the method was validated experimentally. 29 Their results show that a priori knowledge of the location of inhomogeneities is important to know in the reconstruction of optical properties of a heterogeneous tissue. More recently, Sassaroli et al. 79 proposed a fast pMC method for photon migration in a tissue model with an arbitrary distribution of optical properties. This method imposes a minimal requirement on hard disk space; thus it is particularly suitable to solve inverse problems in imaging, such as DOT. Zhu et al. 35 proposed a hybrid approach combining the scaling method and the pMC method to accelerate the MC simulation of diffuse reflectance from a multilayered tissue model with finite-size tumor targets. Besides the advantage in speed, a larger range of probe configurations and tumor models can be simulated by this approach compared to the scaling method or the pMC method alone.

Hybrid MC Methods
Hybrid MC methods incorporate fast analytical calculations such as diffuse approximation into a standard MC simulation. Flock et al. 80 proposed a hybrid method to model light distribution in tissues. In this model, a series of MC simulations for multiple sets of optical properties and geometrical parameters were performed to create a coupling function. Then, this coupling function was used to correct the results computed by diffusion theory. Wang et al. 81 proposed a conceptually different hybrid method to simulate diffuse reflectance from semi-infinite homogeneous media. Wang's method combined the strength of MC modeling in accuracy at locations near the light source and the strength of diffusion theory in speed at locations distant from the source. Wang et al. 82 later extended this method from semiinfinite media to turbid slabs with finite thickness, which is more useful than the previous method in practice. Alexandrakis et al. 62 proposed a fast diffusion-MC method for simulating spatially resolved reflectance and phase delay in a two-layered human skin model, which facilitates the study of frequency-domain optical measurements. This method has been proven to be several hundred times faster than a standard MC simulation.
Hayashi et al. 83 presented a hybrid method to model light propagation in a head model that contains both high-scattering regions and low-scattering regions. Light propagation in high-scattering regions was calculated by diffusion approximation and that in the low-scattering region, i.e., the cerebrospinal fluid layer, was simulated by the MC method. Since the time-consuming MC simulation is employed only in part of the head model, the computation time is significantly shorter than that of the standard MC method. Donner et al. 84 presented a diffusion-MC method for fast calculation of steady-state diffuse reflectance and transmittance from layered tissue models. In their method, the steady-state diffuse reflectance and transmittance profiles of each individual layer were calculated and then convolved to generate the overall diffuse reflectance and transmittance to eliminate the need of considering boundary conditions. Luo et al. 85 introduced an improved diffusion model derived empirically. Then the modified diffusion model was combined with the MC method to estimate diffuse reflectance from turbid media with a high ratio of the absorption coefficient to the reduced scattering coefficient, which can be as large as 0.07. Di Rocco et al. 86 proposed a hybrid method to speed up MC simulations in slab geometries including deep inhomegeneities. In this approach, the tissue model was treated as two sections, i.e., the top layer with a thickness of d in which there is no inhomogeneity and the bottom layer with inhomgeneity. Propagation up to the given depth d, i.e., the top layer, is replaced by analytical calculations using diffusion approximation. Then photon propagation is continued inside the bottom layer using MC rules until the photon is terminated or detected. Tinet et al. 58 adapted the statistical estimator technique used previously in the nuclear engineering field to a fast semi-analytical MC model for simulating time-resolved light scattering problems. There were two steps in this approach. The first step was information generation, in which the contribution to the overall reflectance and transmittance was evaluated for each scattering event. The second step was information processing, in which the results of first step were used to calculate desired results analytically. Chatigny et al. 87 proposed a hybrid method to efficiently model the timeand space-resolved transmittance through a breast tissue model that was divided into multiple isotropic regions and anisotropic regions. In this hybrid method, the standard MC method incorporated with the isotropic diffusion similarity rule was applied to the area that contains both isotropic and anisotropic regions, while the analytical MC, which is similar to Tinet's method, was used for the area that contains isotropic regions only.

Variance Reduction Techniques
In addition to hybrid methods reviewed above, multiple variance reduction techniques, which were initially applied in modeling neutron transport, 88 have also been investigated in the MC modeling of light transport in tissues. For example, the weighted photon model and Russian roulette scheme have been employed in the public-domain MC code, Monte Carlo modeling of photon transport in Multi-Layered tissues. 8 Liu et al. 89 have used one of the oldest and the most widely used variance reduction techniques in MC modeling, i.e., geometry splitting, to speed up the creation of an MC database to estimate the optical properties of a two-layered epithelial tissue model from simulated diffuse reflectance. In this strategy, the tissue model is separated into several volumes, and the technique can reduce variances in certain important volumes by increasing the chance of sampling in important volumes and decreasing the chance of sampling in other volumes. Chen et al. 90 proposed a controlled MC method in which an attractive point with an adjustable attractive factor was introduced to increase the efficiency of trajectory generation by forcing photons to propagate along directions more likely to intersect with the detector, which is similar to geometry splitting in principle. They first demonstrated this approach in transmission geometry 90 and then in reflection geometry. 91 Behin-Ain et al. 92 extended Chen's method for the efficient construction of the early temporal PSF created by the visible or near-infrared photons transmitting through an optically thick scattering medium. More recently, Lima et al. 93,94 incorporated an improved importance sampling method into a standard MC for fast MC simulation of time-domain OCT, by which several hundred times of acceleration has been achieved.

Parallel Computation-Based MC Methods
Parallel computation has received increasing attention recently in the study of speeding up MC simulations due to advances in computer technology. The acceleration due to parallel computation is independent of all previous techniques and thus could be used in combination with them to gain extra benefit. Kirkby et al. 95 reported an approach by which one can run an MC simulation simultaneously on multiple computers, aiming to utilize the unoccupied time slots of networked computers to speed up MC simulations. This method has reduced simulation time appreciably. However, it can be time-consuming to wait for all computers to update the result files in order to get the final result. Moreover, the requirement of saving disk space imposes the use of binary files, and this raised compatibility issues across in various types of computers. Colasanti et al. 96 explored a different approach to address the limitations associated with Kirkby's method. They developed an MC multiple-processor code that can be run on a computer with multiple processors instead of running on many single-processor computers. The results showed that the parallelization reduced computation time significantly.
Considerable efforts have also been made to implement MC codes in graphics processing unit (GPU) environment to speed up MC simulations. Erik et al. 97 proposed a method that was executed on a low-cost GPU to speed up the MC simulation of time-resolved photon propagation in a semi-infinite medium. The results showed that GPU-based MC simulations were 1000 times faster than those performed on a single standard central processing unit (CPU). The same group 98 further proposed an optimization scheme to overcome the performance bottleneck caused by atomic access to harness the full potential of GPU. Martinsen et al. 99 implemented the MC algorithm on an NVIDIA graphics card to model photon transport in turbid media. The GPU-based MC method was found to be 70 times faster than a CPU-based MC method on a 2.67 GHz desktop computer. Fang et al. 100 reported a parallel MC algorithm accelerated by GPU for the simulation of time-resolved photon propagation in an arbitrary 3-D turbid media. It has been demonstrated that GPU-based approach was 300 times faster than the conventional CPU approach when 1792 parallel threads were used. Ren et al. 40 presented an MC algorithm that was implemented into GPU environment to model light transport in a complex heterogeneous tissue model in which the tissue surface was constructed by a number of triangle meshes. The MC algorithm has been tested and validated in a heterogeneous mouse model. Leung et al. 101 proposed a GPU-based MC model to simulate ultrasound modulated light in turbid media. It was found that a GPU-based simulation was 70 times faster compared to CPU-based approach on the same tissue model. Most recently, Cai et al. 102 implemented a fast perturbation MC method proposed by Angelo 79 on GPU. It has been demonstrated that the GPU-based approach was 1000 times faster compared to the conventional CPU-based approach.
Besides using GPU to speed up the MC simulations, some researchers have explored using field-programmable gate arrays (FPGA) to accelerate MC simulations. For example, Lo et al. 103 implemented an MC simulation on a developmental platform with multiple FPGAs. The FPGA-based MC simulation was found to be 80 times faster and 45 times more energy efficient, on average, than the MC simulation executed on a 3 GHz Intel Xeon processor.
Recently, Internet-based parallel computation has gained increasing attention for fast MC modeling of light transport in tissues. Pratx et al. 104 reported a method for performing MC simulation in a massively parallel cloud computing environment based on MapReduce developed by Google. For a cluster size of 240 nodes, an improvement in speed of 1258 times was achieved as compared to the single threaded MC program. Doronin et al. 105 developed a peer-to-peer (P2P) MC code to provide multiuser access for the fast online MC simulation of photon migration in complex turbid media. Their results showed that this P2P-based MC simulation was three times faster than the GPU-based MC simulations.

Acceleration of MC Simulation of Fluorescence
The methods reviewed above are all about the acceleration of MC simulation of diffuse reflectance or transmittance. Compared to diffuse reflectance, fluorescence simulation is more complex and much more time-consuming due to the generation of fluorescence photons upon each absorption event of an excitation photon. A number of groups [11][12][13]30,57,[106][107][108] have employed MC modeling to simulate fluorescence in tissues due to the growing interest in fluorescence spectroscopy or imaging for medical applications. [109][110][111][112] As a consequence, multiple groups have investigated various methods to speed up the MC simulation of fluorescence in biological tissues. Swartling et al. 13 proposed a convolution-based MC method to accelerate the simulation of fluorescence spectra from layered tissues. Their method exploited the symmetry property of the problem, which requires the multilayered tissue model to be infinite in the radial dimension. Different from the conventional fluorescence MC code, this method computed the excitation and emission light profiles separately, from which the spatial distribution of absorption and emission probabilities were obtained. Then a convolution scheme will be applied on the absorption probability and emission probability data to get the final fluorescence signals. Swartling's method has been used by Palmer et al. 113,114 to create an MC database for fluorescence spectroscopy to estimate the fluorescence property of a breast tissue model from fluorescence measurement using a fiber-optic probe. Liebert et al. 57 developed an MC code for fast simulation of time-resolved fluorescence in layered tissues. In this method, both the spatial distribution of fluorescence generation and the distribution of times arrival (DTA) of fluorescence photons at the detectors were calculated along the excitation photons' trajectories. Then the distribution of fluorescence generation inside the medium and DTA as well as the fluorescence conversion probability were used to calculate the final fluorescence signal. It should be noted that the reduced scattering coefficients at the excitation and emission wavelengths have to be approximately equal in this method.

Comparison of Methods for MC Acceleration
Most methods surveyed in the previous sections have been compared and summarized in Table 1 with respect to their acceleration performance, relative error in simulated optical measurements, respective advantages, and limitations. It should be noted that those parallel computation-based methods were not listed in this table because its performance highly depends on the computing architecture, and all the methods summarized in this table can be further accelerated by applying parallel computation.

Applications of MC Methods in Tissue Optics
The most common application of MC method in tissue optics is the simulation of optical measurements such as diffuse reflectance, transmittance, and fluorescence for a given tissue model and illumination/detection geometry, which is considered as a forward problem. In this situation, MC simulations could provide guidelines for the selection of optimal illumination/ detection geometry for selective optical measurements. 46,49,[115][116][117][118] In contrast, MC simulations can also provide data to estimate the optical properties of a tissue model from optical measurements, which is considered as an inverse problem. Solving an inverse problem typically involves the use of a nonlinear least square error algorithm 47,89 or a similar algorithm to find a set of optical properties that would yield optical measurements in MC simulations best matching the actual measurements. Due to the slow speed of traditional MC simulations, a database is frequently created a priori in such an inverse problem to speed up the inversion process. 119 Most of the acceleration methods discussed above can be employed in the creation of such an MC database.
The MC method has been frequently used to find the optimal optical configuration in LDF, one of the oldest techniques in biomedical optics during the past decade. Jentink et al. 120,121 used MC simulations to investigate the relationship between the output of laser Doppler perfusion meters and the optical probe configuration as well as the tissue scattering properties. Stern et al. 122 used MC modeling to simulate the spatial Doppler sensitivity field of a two-fiber velocimeter, by which an optimal fiber configuration was identified. Similar applications can also be found in Refs. 123 through 125. Recently, MC method has been incorporated into LDF to estimate blood flow 126,127 or the phase function of light scattering. 128 The MC method plays an important role in the selection of optimal configuration for PDT because it can generate light distribution in a complex tissue model for PDT dosage determination. Barajas et al. 129 simulated the angular radiance in tissue phantoms and human prostate model to characterize light dosimetry using the MC method. Liu et al. 130 used the MC method to simulate the temporal and spatial distributions of ground-state oxygen, photosensitizer, and singlet oxygen in a skin model for the treatment of human skin cancer. Valentine et al. 131 simulated in vivo protoporphyrin IX (PpIX) fluorescence and singlet oxygen production during PDT for patients with superficial basal cell carcinoma. Later, the same group 132 used the MC method to identify optimal light delivery configuration in PDT on nonmelanoma skin cancer.
The MC method has also been investigated to simulate the OCT signals 133,134 and images [135][136][137] during past years due to its flexibility and high accuracy. Moreover, with the development of efficient MC methods, researchers have started to explore the MC method for image reconstruction in DOT. 138,139 5 Discussion on the Potential Future Directions Due to advances in computing technology, it is expected that the applications of the MC method will be expanded in the near future. A few potential directions in the development of the MC method are discussed below.

Phase Function of Raman Scattering
Raman spectroscopy has been explored extensively for tissue characterization 15,17,140,141 including cancer diagnosis. 14,[142][143][144][145][146][147][148][149] Depending on whether the excitation light is coherent or incoherent, Raman scattering can be broken down into two categories, i.e., spontaneous Raman scattering or coherent Raman scattering. The signal generated out of spontaneous Raman scattering is typically very weak, in which the probability of generating a Raman photon for every excitation photon is lower 150,151 than 10 −7 . Different from that, coherent Raman techniques utilize laser beams at two different frequencies to produce a coherent output, which result in much stronger coherent Raman signals compared to spontaneous Raman scattering. Because of the high chemical specificity of Raman spectroscopy, it is anticipated that there will be more studies using the MC method for Raman spectroscopy to optimize experimental setup. One important issue in these studies is that, the phase function of Raman scattering from biological components in tissues have not been systematically studied. Recent MC studies on Raman scattering 14,15 assumed isotropic Raman emission. This assumption should work fine for spontaneous Raman scattering according to a numerical study. 152 However, this assumption is not valid for coherent Raman scattering since the angular distribution of Raman emission is affected by both the wavelength of the pump light source and the propagating beam geometry. [152][153][154] A systematic study on the phase function of Raman scattering on the molecule level for Raman active biological molecules such as protein and DNA and on the subcellular level for organelles such as mitochondria will be very helpful, in which one or a couple of key parameters similar to the anisotropy factor in elastic scattering could accurately describe the angular distribution of Raman scattering in most common cases. The use of such validated phase functions in MC simulations will yield more useful information than the simplistic treatment in the current literature.

Incorporation of More Realistic Elastic Light Scattering Model into the MC Method
Despite the exploration of various inhomogeneous tissue models discussed above, including the multilayered tissue model, voxelbased and mesh-based tissue models, these tissue models are all based on a few simple optical coefficients including the scattering coefficients and anisotropy factor to characterize optical scatterers. A complete phase function could be used to provide the comprehensive information related to the morphology of optical scatterers, but it is inconvenient for use and its physical meaning is not straightforward. From these scattering properties, the scatterer size and density can be derived 47,155 if they are assumed to be uniformly distributed spheres with homogeneous density. In many scenarios, these assumptions are not valid. For example, it is commonly known that the size and shapes of cells vary significantly with the depth from the tissue surface, and they also change with carcinogenesis. From this point of view, the superposition of multiple phase functions 156 or the fractal distribution of the scatterer size 157 have been proposed to accommodate special situations. An equiphase-sphere approximation for light scattering has also been proposed by Li et al. 158 to model inhomogeneous microparticles with complex interior structures. Later, the same group reported two stochastic models, 159 i.e., the Gaussian random sphere model and the Gaussian random field model, to simulate irregular shapes and internal structures in tissues. The incorporation of these more realistic elastic light scattering models into the MC method will expand its capability and offer more accurate information about light scatterers in tissues.

Exploration of the MC Method in Imaging Reconstruction
In most current applications of the MC method, the tissue model is assumed to be a simple layered model or determined a priori, which does not fully exploit the potential of the MC method in preclinical or clinical imaging/spectroscopy. When the MC method becomes adequately fast in the near future, which might be mostly attributed to the combination of the accelerated MC methods and parallel computing discussed above, the MC method could be used to reconstruct the optical properties of a complex tissue in optical tomography, in which the morphological structure could be obtained in real time by fast imaging techniques such as OCT for superficial regions of interest or magnetic resonance imaging for deep regions of interest in a large tissue volume. Those voxelated tissue models, 19,[36][37][38][39][40][41][42] which are more realistic than simplified homogeneous or layered tissue models, can be readily used in such reconstruction. The speed of reconstruction might be comparable to that using diffusion approximation reported in the current literature, but the accuracy would be considerably higher in a tissue model on the millimeter scale.

Conclusion
In this review, the principles of MC modeling for the simulation of light transport in tissues were described at the beginning. Then a variety of methods for speeding up MC simulations were discussed to overcome the time-consuming weakness of MC modeling. Then the applications of MC methods in biomedical optics were briefly surveyed. Finally, the potential directions for the future development of the MC method in tissue optics were discussed. We hope this review has achieved its intended purpose to give a general survey on the capability of MC modeling in tissue optics and other relevant key techniques for those readers who are interested in MC modeling of light transport.