Parallelized Monte Carlo Photon Transport Simulations for Arbitrary Multi-Angle Wide-Field Illumination in Optoacoustic Imaging

Optoacoustic imaging (OAI) or photoacoustic imaging can resolve the distribution of tissue chromophores and optical contrast agents deep inside tissue from the optoacoustic detection with multi-spectral illumination. The development of a fast and accurate modeling method for the photon transport in OAI is necessary to quantitatively evaluate the tissue optical parameters. This paper presents a parallelized Monte Carlo modeling method especially for OAI (MCOAI) to simulate photon transport in bio-tissues with arbitrary multi-angle wide-field illumination. The performance of the MCOAI method is verified by comparison with the graphics processing unit (GPU)-accelerated MCX method in the typical cases with the pencil beam and ring source illumination. The simulation results demonstrate the GPU-based MCOAI method has equivalent accuracy and significantly improved computation efficiency, compared to the MCX method. The simulations with cylindrical and hemispherical source illumination further illustrate that the MCOAI method can effectively implement three-dimensional photon transport simulation for typical illumination geometries of OAI systems. A cross-section of Digimouse is selected as a realistic heterogeneous phantom illuminated with six different light sources usually employed in OAI, in order to prove the necessity of establishing photon transport modeling in OAI for the quantitative visualization in deep tissues. The MCOAI method can provide a powerful tool to efficiently establish photon transport modeling with arbitrary illumination modes for OAI applications.


INTRODUCTION
Optoacoustic imaging (OAI) or photoacoustic imaging is a hybrid imaging modality that can provide multi-scale optical visualization with high depth-to-resolution ratio in biological tissues [1][2][3][4][5][6]. The distribution of the initial acoustic pressure or absorbed light energy density P 0 , is generally reconstructed from optoacoustic signals detected by ultrasonic transducers in traditional OAI methods. The initial acoustic pressure P 0 is the product of optical absorption coefficient and light fluence in the irradiated region [7,8]. Since light fluence in the tissue is attenuated along the penetration depth, the result of traditional OAI, i.e., the distribution of P 0 , cannot directly reflect the optical properties of deep tissues, such as optical absorption coefficient [9,10]. To accurately reveal the optical functional information of deep tissues, a precise photon transport modeling approach is necessary to be developed for offering assistance to quantitative OAI methods.
As the gold standard for predicting distribution of light fluence in the biological tissue, the Monte Carlo (MC) simulation with strong applicability is preferred for the precise modeling with multi-angle wide-field illuminations usually employed in OAI systems. MC simulation methods have been widely developed for optical imaging modalities [11][12][13][14][15][16]. The MCML and MCXYZ methods by Wang et al. have been used to study light propagation in multi-layered tissues [11,12]. Several light source modes have been realized including the pencil beam and the finite beam with a fixed incident direction, as well as the isotropic beam. An accelerated MCML has been developed to reduce the computation time based on graphics processing unit (GPU) [13]. The MMC and GPU-based MCX methods by Fang et al. provide more light source modes, containing the slit beam, the planar beam, and the cone beam, with the option of sinusoidal modulated amplitude [14,15]. For the multiangle wide-field illumination modes usually employed in OAI systems, the conventional MC methods can not simultaneously and randomly generate the initial angles and positions of emitted photons using the GPU framework. Thus, the number of initial angles and positions is limited due to the low computation efficiency of the central processing unit (CPU) framework. The limited number of initial positions produces great challenges for efficient light transport modeling of arbitrary multi-angle widefield illumination modes in OAI systems, such as ring, cylindrical and hemispherical sources for typical OAI systems.
In this work, we present a MC approach based on GPU to quickly and accurately establish light transport modeling for OAI systems with arbitrary multi-angle wide-field illuminations, termed as MCOAI. The incident angles and positions of emitted photons can be simultaneously and randomly generated from their pre-set probability distributions in MCOAI, allowing its effective performance for simulating arbitrary light sources with GPU acceleration. By comparing with the MCX method, the accuracy of MCOAI method is verified in the simulation with the pencil beam illumination. Comparison results using the ring source mode in optoacoustic tomography indicate that the MCOAI method performs higher efficiency and greater precision than the MCX method. The proposed MCOAI method can also effectively achieve three-dimensional (3D) light transport simulation for cylindrical and hemispherical illumination modes usually employed in OAI systems. Finally, several illumination modes are set up for a cross-section of Digimouse, illustrating the necessity of MCOAI for the quantitative OAI.

MCOAI Method
The MCOAI method can achieve simultaneous generation of both initial angles and positions of emitted photons by random numbers from pre-set probability distributions, or random selection from customized tables in the GPU framework.
Owing to the sufficient number of initial angles and positions generated in the GPU framework, the MCOAI method can effectively conduct the precise photon transport modeling of arbitrary 3D media illuminated by arbitrary multi-angle widefield light sources. In Cartesian coordinates, all objectives are uniformly divided into voxels. Optical properties of the media, such as absorption coefficient µ a , scattering coefficient µ s , scattering anisotropy factor g and refractive index n, are assigned to each voxel. The incident beam is simulated by emitting photon packets, whose trajectories are calculated in the media. Probability models and pseudorandom number generators are utilized to determine the initial position and direction, the path length between two adjacent scattering events, the scattering angle, transmission or reflection, as well as termination of the photon packet. As photons propagate from one voxel to the next, they deposit some energy ("weight") into the voxel, depending on absorption coefficient µ a of the voxel. Photons are tracked until terminated in Russian roulette or escaping from the media boundary. The deposited energy is digitally accumulated in a 3D matrix, and is normalized to output power at the end of the simulation. Figure 1 illustrates the flowchart of MCOAI for the multi-angle wide-field illumination mode.
In MCOAI, the random number is randomly generated in [0, 1]. For the ring source with the radius of r in the plane of z = z 0 , a random number R 0 is used to obtain the longitude angle φ = 2πR 0 uniformly distributed in [0, 2π]. After an initialized direction φ is randomly generated, the initialized position r s and the direction vector − → v s can be expressed as For the cylindrical source with the radius of r and the height of H, a random number R 1 is employed to produce the longitude angle φ = 2πR 1 , and another random number R 2 is used to compute the elevation position h = HR 2 randomly distributed in [0, H]. The initialized position r s and the direction vector − → v s can be expressed as For the hemispherical source with the radius of r, two random numbers R 3 and R 4 are used to generate the longitude angle φ = 2πR 3 and latitude angle θ = πR 4 /2, respectively. The initialized position r s and direction vector − → v s can be expressed as For other types of source mode, the initial direction and position of emitted photon can be customized aforehand, and then the random extraction is adopted to randomly obtain these parameters based on the probability distributions. For light source models with uneven light intensity distribution, packet weights at different

Simulation Set-Up
A total of 10 8 photons are launched in the simulation of the proposed MCOAI method. One hundred and twenty-eight threads are set for parallel calculation. Four simulated samples are employed to validate the performance of the MCOAI method. The optical properties of the samples are set as the representative values in the optical window [16]. The MCOAI and MCX methods are run on the same computer with an NVIDIA GTX 750 GPU with 64 GB memory.

1) Sample 1: A homogeneous cube:
A cubic domain of 60 mm × 60 mm × 60 mm is simulated shown in Figure 2. The volume is from −3 to +3 cm in x and y axis, and from 0 to 6 cm in z axis, with the number of grids in each axis of 60. The pencil beam is located at r s = (30,30,0) with an incident direction vector − → v s = (0,0,1); the medium is homogeneous with the absorption coefficient µ a = 0.01 mm −1 , the scattering coefficient µ s = 10 mm −1 , the refractive index n = 1.37 and the anisotropy g = 0.9.
2) Sample 2: A homogeneous cylinder: A cylinder sample with the center located at (6,6,8) mm has a radius of 5 mm and a height of 6 mm. The voxel size is 50 µm × 50 µm × 50 µm to build the accurate curved boundary. The optical properties for the cylinder are µ a = 0.1 mm −1 , µ s = 10 mm −1 , n = 1.37 and g = 0.9. A ring source with the radius of 5 mm is located at the plane of z = 8 mm, shown in Figure 3A. Figure 4A illustrates a 3D cylindrical illumination distributed around the cylinder boundary.
3) Sample 3: A homogeneous hemisphere: A hemisphere sample centered at (6, 6, 6) mm has a radius of 5 mm. The voxel size and optical properties of the hemisphere are the same as those of the Sample 2. Figure 4C shows a hemispherical illumination distributed around the hemisphere surface. 4) Sample 4: Cross-section of Digimouse: A cross-section of the widely-used Digimouse atlas has been chosen as a more realistic tissue sample. g = 0.9 and n = 1.37 are set to be homogeneous in the sample. µ a and µ s of the sample are shown in Table 1.
Several types of light sources are simulated, including the finitewidth light-sheet illumination from one side (left, right, top, or bottom) and four sides, as well as ring light-sheet illumination.

RESULTS AND DISCUSSION
To evaluate the accuracy of the MCOAI method, a comparison with the MCX method is implemented using the pencil beam to illuminate Sample 1. The fluence maps at the plane of y = 30 mm from the MCOAI and MCX results are presented in Figures 2A,B, respectively. The fluence contours with 10 dB spacing are also displayed in Figure 2C. The MCOAI result matches well with the MCX result in both fluence maps and  contours, presenting the accurate performance of the MCOAI method. To quantitatively assess the accuracy of MCOAI result in Figure 2, the root mean square error (RMSE) has been calculated to show the difference between the MCOAI result and the MCX result. The RMSE is defined as where N is the total number of the image pixels, r MCX (n) and r MCOAI (n) are the values of fluence maps at the nth pixel in the MCX result and MCOAI result, respectively. The RMSE value between Figures 2A,B is 0.0008, indicating high consistency of the MCOAI result and the MCX result.
To demonstrate the fast and accurate modeling of the MCOAI method for the multi-angle wide-field illumination mode, a cylindrical sample (Sample 2) in Figure 3A has been chosen, with a ring source illumination usually employed in OAI. The calculation time of the MCOAI and MCX methods is illustrated in Table 2. The run time of MCOAI is 1.7 min. In the MCX method, amounts of pencil beam sources are set at the assigned positions around the cylinder boundary. The MCX simulations     Figure 3E is obviously smoother and more authentic than the images in Figures 3F,G. Since 10 8 photon packets are emitted at random initialized positions along the ring shape in the MCOAI method, the finite number of initialized positions set in the MCX method is much less than the number in the MCOAI method, resulting in the uneven fluence distributions at the boundary in the MCX results. The comparison results of the focused light illumination in Figure S1 are consistent with the results of the ring source in Figure 3. The standard deviation (SD) has been calculated to analyze the regional uniformity of the MCOAI and MCX results indicated by the black squares in Figures 3E-G. The SD is defined as where r(n) is the value of fluence map at the nth pixel and r mean is the mean value in the black square. The SD values of the MCOAI result, the MCX result of 360 sources and the MCX result of 720 sources are 0.0142, 0.0828, and 0.0520, respectively. The lowest SD value of the MCOAI result indicates the highest uniformity. The higher fidelity and efficiency of the MCOAI method can be attributed to introducing the random generation approach for the incident angle and position of each emitted photon into the parallel computing framework. Therefore, MCOAI method is more appropriate to accurately and quickly establish multi-angle wide-field illumination modes for optoacoustic tomography. Furthermore, two typical 3D source modes (Figures 4A,E) usually existed in OAI are illustrated by MCOAI method. For the cylindrical illumination mode, the fluence maps at the planes of z = 9 mm, y = 6 mm and 3D view are presented in Figures 4B-D. For the hemispherical illumination mode, the fluence maps at the planes of z = 6 mm, y = 6 mm and 3D view are displayed in Figures 4F-H. The boundary and interior distributions of fluence maps in two homogeneous media are smooth and authentic owing to the sufficient incident positions generated randomly. The results demonstrate the light transport modeling in 3D simulation with multi-angle wide-field illumination modes can be efficiently and accurately established by the proposed MCOAI method.
This MCOAI approach is further used to simulate various light illumination modes on a cross-section of Digimouse (Sample 4) shown in Figure 5. The optoacoustic tomography technology has been widely employed for the whole-body smallanimal imaging by integrating lots of cross-section images. The cross-section of Digimouse containing several main organs listed in Table 1 is selected.
The MCOAI method was further utilized to illustrate the necessity for quantitatively evaluating functional information of biological tissue from OAI results. Figures 5C,D present the grayscale images of µ a and µ s of Sample 4, indicating the heterogeneity in different organs. Figure 5E shows the images of initial acoustic pressure P of Sample 4 using six different illumination modes. As displayed in the first four columns of Figure 5E using one-sided illumination, the intensity distribution of targets and background in the P images plummet when moving away the light source. In the last two columns of Figure 5E, due to the illumination modes are both full-angle wide-field, the P images are closer to the µ a images while still presenting obvious differences. The results show the P image cannot directly reflect accurate quantitative distributions of µ a in the sample, affected by the light attenuation along tissue depth.
The fluence maps under these six different illumination modes are shown in Figure 5F, and the contours with 10 dB spacing are also displayed in Figure 5G. The light fluence φ in deep tissues is attenuated along the penetration depth, thus the fluence distributions in the region near the light source are higher than those in the region far from the light source. Under the finitewidth light-sheet illumination modes only from one side, the heterogeneous distribution of φ is extremely evident. Under the four-sided light-sheet and ring light-sheet illumination modes, the Sample 4 receives full-angle wide-field irradiation, the fluence maps φ are attenuated from boundary to center. The results demonstrate the MCOAI method should be utilized to provide the fluence distribution of φ, in order to quantitatively recover the distribution of optical properties.

CONCLUSION
To quantitatively evaluate the optical properties of the biological tissues, it becomes essential to develop an efficient and precise photon transport modeling method for multi-angle wide-field illumination modes usually employed in OAI systems. In this work, we present the MCOAI method to randomly generate the incident angles and positions of emitted photons based on parallel computing framework for efficiently simulating arbitrary wide-field light sources. The MCOAI method can offer the same accuracy in the simulation using pencil beam as the MCX method. However, the MCOAI method for the multi-angle wide-field light sources, such as the ring source, is capable of providing more accurate and authentic distributions of light fluence in a few tenths of the time, compared to the MCX method. The GPU framework has been utilized in both the MCX method and the MCOAI method for acceleration. For the pencil beam and planar beam in the MCX method, the initial angles and positions of emitted photons are generated in the GPU framework. However, for the multi-angle widefield illumination modes usually employed in OAI systems, the MCX method can not simultaneously and randomly generate the initial angles and positions of emitted photons in the GPU framework. The proposed MCOAI method can achieve simultaneous generation of both initial angles and positions of emitted photons by random numbers in the GPU framework, greatly reducing the running time. If the optical properties of samples at the selected wavelength in visible or near-infrared spectrum can be determined, the light transport modeling at any selected wavelength can be achieved by the proposed MCOAI method. The proposed method has been demonstrated with typical illumination modes and uniform light intensity, which can be explored on user-defined illumination modes and uneven light intensity for further study. The MCOAI method is able to precisely conduct the light transport modeling for all typical OAI systems, having great potential of contributing to the further investigation of quantitative OAI methods. In this work, the MCOAI method focuses on providing an accurate light transport simulation for the forward modeling in quantitative OAI, not involving the inverse problems. The MCOAI method is designed to be an open source for user-friendly employment.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. The data of this article will be made available by the authors to any qualified researcher.