Medical optical applications require knowledge about the optical properties in target tissue. The radiance transport equation (RTE) is an accurate model for light propagation in tissue media. Monte Carlo simulations (MCs) and Diffusion equations (DE) are two widely used calculation methods to solve the RTE [1-3]. MCS displays advantages over other methods due to its flexibility and accuracy. In MCS, the light dose is computed by simulating a large number of photons’ propagations with random numbers. Done this way, MCS is a time-consuming method. For instance, CPU-based MCS may take hours to complete one accurate calculation.

Recent improvements in the programmability of graphics processing units (GPU) have enabled use of the GPU to accelerate the MC simulations [4-7]. Alerstam et al. presented a GPU-based MC to simulate a homogeneous medium, providing a massive 1000 × speed-up over traditional CPU code [4]. Additionally, Fang’s MC model achieved a 100 × speed-up by using the GPU to trace photon migration in a heterogeneous medium [5]. Moreover, Jiang et al. demonstrated a GPU cluster-based Monte Carlo simulation [8], whose speed up is proportional to the number of GPUs used. Recently, Doronin et al. developed an online object oriented Monte Carlo computational tool [9], which makes the GPU-based simulations more easily accessible for the scientific community. Therefore, based on the rapid development of GPU and internet technology, MCS is becoming a real-time and accessible method.

It is widely acknowledged that MCS is a statistical method that relies on an abundance of repeated random sampling to reduce the statistical error. When the calculation time of MCS is several hours, it is acceptable to take several minutes to determine the number of simulated photons, which is known as batch method [10] to achieve a particular accuracy. However, as the calculation speed issue can be overcome by parallel processing based GPU calculation, the batch method emerges as an inconvenience to utilize MCS. Therefore, a dynamic accuracy criterion, which is also known as history by history statistical estimators [10, 11], is highly demanded. In the next section, a dynamic accuracy estimation for GPU-based Monte Carlo simulation will be described in detail. The simulations are performed on a notebook computer (NVidia GeForce 840M, 8G RAM and 2.0G CPU).

A tissue optical measurement, consists of three basic elements: laser sources, tissue medium and optical sensors [1]. Usually, a laser source, with P₀ mw intensity is simulated by a pencil beam with N₀ photon package in MCS. Each photon packet is initially assigned a weight W, equal to unity. The position and detected surface of the optical sensor are preset according to the real set-up. The tissue medium is described by the following parameters: the absorption coefficient μ_a [cm⁻¹], the scattering coefficient μ_s [cm⁻¹] and the anisotropy factor g. The code and detailed description of the MCS process can be found elsewhere [1-2, 4-5].

Assumptions in MC simulation, after a large number of photon tracing, the sensor detects N photons, and then the probe light intensity value is

denotes the optical weight for the i^th photon detected by the sensor. N is defined as the detected photon number. Based on the Law of Large Numbers, P_det gets closer to theoretical value P_{det 0} as N gets larger. P₀ and N₀ are pre-defined, thus it is critical to derive the statistical error for .

In MCS, if is regarded as a random variable x, whose standard deviation is σ{x} , then the standard deviation for is . As the simulated photon’s number N is large enough, σ²{x} can be approximated by:

Thus, the relative error can be estimated by the following equation:

In a GPU-based MCS, there are Nt threads running parallel simultaneously, then Eq. (3) can be approximated by:

represents the optical weight for the i^th photon detected by the sensor in thread j, N is defined as the detected photon number in thread j.

This work utilizes a widely accepted GPU-based MCS code from Ref. [5]. It is worth mentioning that, Eq. (4) can easily be applied as the dynamic accuracy estimation in the GPU-based MCS code to predict the statistical error dynamically. In the next section, MCS for homogenous tissue medium and heterogeneous (brain) tissue medium will be performed and the dynamic accuracy estimation will be tested in these simulation examples.

To verify the feasibility of the dynamic accuracy estimation, an implementation detail of MCS for a homogenous tissue medium is presented here. The statistical errors, which are predicted by dynamic accuracy estimation and batch method, are also analyzed in this section.

In Fig. 1, a fiber bundle probe is placed on a semi-infinite homogenous tissue medium. In the optical fiber, the light source is incident to the biological medium through fiber F0, which is denoted as source F0, and the optical fibers F1-F4, which collect the back scattered light, are denoted as sensors F1-F4. The core radius of the fiber is 200 μm. The distance of the optical fiber is 280 μm. Hence, the distance between source fiber F0 and detect fiber F4 is 1120 (4×280) μm. The refractive index for the cladding and core are 1.4413 and 1.458, respectively. The absorption coefficient, the scattering coefficient and the anisotropy factor g for the tissue medium are 0.1 cm^-1, 60 cm^-1 and 0.8, respectively.

The distance between sensor F4 and source F0 is the largest; hence F4 suffers the most serious error. When the sensor F4 has a great computational accuracy, it can be considered that the sensors F1-F3 also achieve the ideal calculation precision. Therefore, only the statistical error of sensor F4 is considered herein.

In each thread, there is a global initially assigned variable rep_j, equal to zero. When F4 detected a photon package with weight w, rep_j is modified by the following equation:rep_j = rep_j + (w − w')². w' is the photon package previous detected by F4. After the above process, w' = w and w_{det j} = w_{det j} + w'. Herein, w_{det j} is another global variable, which is utilized to record the total detected photon package’s weight. Subsequently, the relative statistical error can be predicted by the following equation:

For ease of exposition, the statistical error, which is predicted by the dynamic accuracy estimation, denoted as dynamic statistical error.

We also applied the traditional statistical error predicted method in this work. As the sensor F4’s detected photon number equals to a preset number, the MCS process ceases and the results are stored in the hard disk. is denoted as the total detected photon package’s weight by sensor F4 in one simulation. Ten MCS’s are implemented, afterwards, the average detected weight can be calculated by . The statistical error predicted by the batch method [10], which is denoted as static statistical error, can be described by the following equation:

n is equal to 10 and represents ten simulations in this study. When the sensor F4’s detected photon number reaches a preset quantity, the simulation terminates, and the dynamic and the static relative statistical errors are calculated based on the above discussion. The simulation results are shown in Fig. 1(b). Based on the nature of MCS, the error is reduced gradually as the detected photon number increases. Indeed, our results are consistent with this theory. There is a deviation between the dynamic and the static relative statistical error estimation, which is about 0.1% to 0.4%. The deviation stems from the approximation in Eq. (2). In fact, the power stability of a normal laser is about 1%. Therefore, 0.1% to 0.4% deviation exerts minimal effect on the actual application. It’s worth mentioning that the utilization of dynamic accuracy estimation has little influence on the efficiency of an MCS program. The computation time is slightly increased by about 1/69 compared to the traditional MCS without dynamic accuracy estimation.

Light dose computation in the tissue medium is important for photodynamic therapy (PDT) [12] and functional near-infrared spectroscopy (fNIRS) [13]. To improve the efficiency of light dose computation, numerous previous studies had introduced GPU-based MCS in a biomedical model, such as skin [4], brain [5] and mouse [6] numerical optical model. With more computing cores, the GPU computing performance has been steadily increasing and achieves a billion floating-point operations per second (GFLOPS). It is apparent that the MCS will be applied to the clinical field in the near future. In our work, this dynamic accuracy estimation is applied in the MCS for brain numerical model to quickly monitor statistical error and halt the simulation as the statistical error reaches less than 1%.

The brain numerical optical model has been described in detail in Ref. [2] and [5]. In our work, the light is input from the forehead (as indicated by a solid arrow in Fig. 2(a)) and the target tissue is prefrontal cortex, which relates to higher cognitive function, such as memory and decision making. To utilize dynamic accuracy estimation in these simulations, a virtual sensor is placed at white matter (depth of 21 mm, 32 mm apart from the input position, as indicated by a dashed arrow in Fig. 2(a)). Clinically, the virtual sensor can be placed on the region of interest. The logarithm results of light doses inside the brain are shown in Fig. 2. The simulated photon numbers are 10⁵, 10⁷, and 10⁸ for Figs. 2(a)-2(c), and the simulation time is 1.92 s, 164.73 s and 1515.93 s. When the simulated photon number is 10⁵, the calculation accuracy is unsatisfied (Error=4.7%) and the dose distribution image suffers from statistical noise, especially in the vicinity of the virtual sensor. As the simulated photon number increases, the light dose results get more accurate. In Fig. 2(c), the statistical error is reduced to 0.3% and the noise of dose distribution image is also suppressed. Therefore, our results suggest the dynamic accuracy estimation has the ability to real-time monitor the noise in the MCs for light dose calculation. In practices, the simulation can automatically cease as statistical error less than 1%.

Notably, the light is able to penetrate the skull and arrive at the white matter of the brain; hence, the light can be employed to study the surface of the brain in clinical and experimental purposes.

Time-of-flight measurement is another important tool to study blood oxygenation level-dependent (BOLD) signal. In this kind of measurement setup, a pulse laser, with pulse width less than 1 Ps, incident into the head. And a high-speed optical sensor, which is placed next to the laser incident position, monitors the reflectance from the prefrontal cortex. Using the reflectance, one can calculate the blood oxygen concentration of the position in the vicinity of the optical sensor. This technique is known as fNIRS. A time-domain MCS is utilized in our work to simulate a pulse laser propagating in the brain and detected in reflection mode. A virtual sensor is placed at the 1200 Ps time-gate, as indicated by the arrow in Fig. 3(b). When the detected photon number by the sensor is 10⁷, the statistical error is 9.2% and the reflectance curve suffers serious noise (see the sub-figure in Fig. 3(b)). As the detected photon number is increased to 10⁸ or 10⁹, the statistical error is reduced to 2.8% or 0.8% respectively. The reflectance curve is very smooth as shown in Fig. 3(d), which illustrates that the statistical noise is basically cleared. Therefore, the dynamic accuracy estimation can also be applied in time-domain MCS to precisely monitor the noise. It should be noted that the simulation times for results in Figs. 3(b)-3(d) are 135.76 s, 1352.21 s and 13521.57 s.

In our work, a method named dynamic accuracy estimation has been developed to monitor statistical error in the GPU-based MCS. This method is very compatible with both the spatial and time-domain MCS. By using the dynamic accuracy estimation, the computation time of MCS is slightly increased by about 1/69 compared to the traditional one. Through using GPU to accelerate the MC simulations, we can acquire simulation results (statistical error<1%) within 13521 s for a time-of-flight measurement. It is worth noting that our GPU-based MCS is still not yet optimized in terms of simulation speed and we expect significant simulation speed-up improvements by using a multi-GPU platform.

As the GPU-based MCS is becoming more effective, the utilization of dynamic accuracy estimation in MCS can reflect the current calculation accuracy in real time. In addition, the skin and brain are both heterogeneous media, but with different optical properties. If we utilized the optical properties of skin in the program, the method can also be applied for skin simulation. Therefore, our method can not only be used in fNIRS, but also in photodynamics therapy for skin disease.