MCOCT: an experimentally and numerically validated, open-source Monte Carlo simulator for optical coherence tomography

Here, we present MCOCT, a Monte Carlo simulator for optical coherence tomography (OCT), incorporating a Gaussian illumination scheme and bias to increase backscattered event collection. MCOCT optical fluence was numerically compared and validated to an established simulator (MCX) and showed concordance at the focus while diverging slightly with distance to it. MCOCT OCT signals were experimentally compared and validated to OCT signals acquired in tissue-mimicking phantoms with known optical properties and showed a similar attenuation pattern with increasing depth while diverging beyond 1.5 mm and proximal to layer interfaces. MCOCT may help in the design of OCT systems for a wide range of applications.


Introduction
Optical coherence tomography (OCT) is a high-resolution, non-invasive imaging technique increasingly used for medical and biological applications.It is based on the principle of low-coherence interferometry, which relies on utilizing a light source with a broad wavelength bandwidth to ensure coherence gating, thus leading to a high axial resolution of a few microns [1].In OCT systems, an interference spectrum is created by combining the reflected light from a reference mirror with that from the sample arm.This pattern is used to obtain a reflectance profile of the sample as a function of depth, namely an OCT A-line.OCT imaging primarily relies on detecting ballistic and quasi-ballistic backscattered photons from turbid environments, such as biological tissues.As light penetrates deeper into the tissue, its intensity decreases exponentially due to scattering and absorption, subsequently limiting the imaging depth to a few millimeters.OCT is now a standard diagnostic tool in the field of ophthalmology and a promising candidate in various other fields, including cardiology, dermatology, gastroenterology, and laryngology [2][3][4][5].
Light propagation in a turbid medium can be modeled in simple geometries by analytically solving the radiative transfer equation (RTE) [6].However, when imaging biological tissue with high scattering properties, the mean free path of light within the tissue is short, which renders the RTE calculation and its inverse problem reconstruction too complex.In such cases, numerical simulation techniques using Monte Carlo (MC) methods have emerged as the gold standard for modeling photon interactions in biological tissue [7][8][9][10][11][12].With MC methods, light is considered a collection of photon packets tracked as they undergo multiple events of reflection, refraction, absorption, and scattering while traveling through the turbid medium.Over the past years, several MC algorithms have been developed, amongst which Monte Carlo multi-layered (MCML) is the most widely known, forming the foundation of numerous other simulators for biomedical applications [9].In MCML, an ensemble of photon packets is launched as a pencil beam and perpendicularly to the surface of a multi-layered medium, with a simulated beam geometry that is cylindrically symmetric around the z-axis.The use of voxels, tetrahedrons, and 3D volumetric meshes further allowed the extension of planar and cylindrical geometries to model more realistic tissue structures [11,[13][14][15][16].To accelerate the speed and improve computational efficiency, it was also proposed to use a graphics processing unit (GPU), enabling simultaneous simulations of a higher number of photon packets [17,18].
Monte Carlo methods were previously applied to simulate OCT signals [10,[19][20][21][22][23][24][25][26] and some open-source OCT simulators recently became available [27][28][29].Monte Carlo simulation for tissue layers with embedded objects for OCT (MCEO-OCT) is a modified version of MCML that embeds geometric shapes in the multi-layered tissue to mimic vessels and utilizes the method of importance sampling to accelerate the computation time [27].The OCT massively parallel simulator (OCT-MPS) is based on an MC approach that exploits a tetrahedron mesh and GPU programming [28].These two simulators were developed for time-domain OCT applications as only photons describing a path length difference with the reference path length smaller than the coherence length of the OCT system are tagged as contributing to the signal.For Fourier-domain applications, three MC simulators have previously been developed: one uses a layered tissue geometry with Gaussian beam focusing [30], another exploits a 3D voxel system with a pencil beam [31], while the third tool proposed an acceleration scheme based on importance sampling with a phase-function that depends on the wavelength [29].The latter tool is open-source.These tools can be used to understand OCT signals better, aid in the instrumental design and development of OCT-based applications, guide the development of OCT signal processing algorithms, and more recently train neural networks to improve OCT image reconstruction [10,[19][20][21][22][23][24][25]27,28,[30][31][32][33].
In this work, we propose an open-source Fourier-domain simulator, MCOCT, modified from MCML 3D voxel geometry mcxyz, in which we implemented modifications to accommodate the particular requirements of Fourier-domain OCT.To achieve depth-resolved reflectometry in addition to simulating photon propagation, photon packets are simulated at a single central wavelength and tracked to record their traveled distance and extract phase as a function of wavelength.This strategy allows mimicking coherence gating.A photon detection system with an angle of acceptance was also implemented to mimic the numerical aperture of the imaging system [10].To achieve a significant signal-to-noise ratio and expedite computation time, importance sampling was also employed [24,34].Gaussian illumination was further incorporated to model the light intensity profile and better represent OCT imaging.For numerical validation, the optical fluence rate generated by MCOCT was compared to the one simulated by Monte Carlo eXtreme (MCX) [18] using the same optical properties and geometry.For experimental validation, OCT signals generated by MCOCT were compared to OCT signals acquired at specific numerical apertures in synthetic phantoms with known optical properties.

Material and methods
In the following subsections, we will sequentially describe the specific parts of the MCOCT simulator.

Monte Carlo photon propagation
A standard MC simulation starts with the launch of a photon packet into a medium that is defined by its index of refraction n, absorption coefficient µ a , scattering coefficient µ s , and anisotropy coefficient g [9].For voxel-based approaches, each voxel is associated with a tissue type and its optical properties.Photon packets randomly propagating through the geometry are associated upon launch with a weight W = 1 representing the packet intensity that decreases as it undergoes events.At the site of an absorption interaction, W is updated according to W = W • exp −(µ a • s), where s is the distance traveled by the photon in that step [9].In the case of a scattering interaction, the diffusion direction is calculated from the azimuthal angle ϕ and the deflection angle θ s [6].ϕ is uniformly distributed in the range [0, 2π], whereas θ s is typically modeled by the Henyey-Greenstein (HG) function f HG that depends on the optical properties of the medium [35,36].Since most biological tissues have an anisotropy coefficient g ∈ [0.90, 0.95] at the clinically applied wavelengths of OCT, the photons are more prone to scatter in the forward direction [37].When W is below a predetermined critical threshold, a roulette process begins in which the photon packet is associated with a probability to continue its propagation, otherwise, the sequence is terminated for this particular photon packet, and W is zeroed [9].The simulation ends upon the termination of all photon packets.

OCT signal simulation
In order to simulate OCT A-lines, we introduced a new beam-focusing method based on Gaussian beams, used statistical methods to increase the probability of photon detection, tracked photons path lengths, and incorporated a system mimicking photon illumination and detection.

Gaussian beam description
In OCT, beam focussing impacts both signal intensity and resolution.The diffraction-limited focal spot of a Gaussian beam at a certain wavelength λ is determined by its beam waist w 0 .Employing the paraxial approximation, the beam waist can be described as w 0 = 2λf /(π • D), relative to the focal length of the imaging lens f , the beam center wavelength λ and the beam diameter at the lens D [38].In OCT imaging, light propagates through several layers of different indices of refraction n i with each layer affecting the axial position and the width of the beam waist, i.e., the beam curvature.Figure 1 illustrates the curved trajectory of such a beam depicting the individual parameters where z = 0 corresponds to the sample surface, described by the red line.The position of the focus in the tissue is given by z f and the distance between the tissue surface and imaging lens by z f ,I , presenting a more general approach compared to assuming the focusing lens to be at z = 0 [38].
To implement this beam curvature in the MC method, the beam was focused by propagating photons along the trajectory by a small step unit s min to the next position (x, y, z).Thereby, the unit of the step size is defined by the optical properties of the tissue.After each step, the direction unit vector û is updated to ensure precise photon tracking along the initial trajectory.Considering a photon propagation through two layers of differing refractive indices n 1 and n 2 , the beam waist of the initial beam in a medium with n 2 can be described as The Rayleigh range of the beam z R in the medium of n 2 is and the radius of curvature of the beam R(z) in the medium at a position z can be expressed as where z f ,L is the axial distance between the focal plane in the medium and the photon position, which is updated after the photon moves.The resulting direction vector has an orientation The photon packet propagates by  min until it undergoes its first scattering event, after which it propagates in accordance with a traditional MC simulation [9].To allow the photon packet to follow a curved path, the packet is moved by a distance  min before its direction û is recalculated using Eq. 4. The photon is moved this way until it has traveled the entire distance .To assess the performance of the Gaussian illumination, simulations were performed with a pencil beam and a geometric focused beam [9,18].Variables are described in Table 5 of the Appendix.

Importance sampling and photon splitting methods
To increase the probability of backscattering events, the method of importance sampling was implemented and combined with photon splitting (see flowchart in Fig. 2) [24,34].Importance sampling is a probabilistic method to accelerate MC simulations for OCT by biasing random events leading to backscattering, thereby increasing the probability of photon detection and thus enhancing computational efficiency.An updated biased scattering function   based on the HG probability density function   was then introduced.After a scattering event, this function modified the photon propagation direction to be oriented towards the detector, i.e., the collecting perpendicular to the radius of the curvature of the beam and is consequently given by The photon packet propagates by s min until it undergoes its first scattering event, after which it propagates in accordance with a traditional MC simulation [9].To allow the photon packet to follow a curved path, the packet is moved by a distance s min before its direction û is recalculated using Eq. ( 4).The photon is moved this way until it has traveled the entire distance s.To assess the performance of the Gaussian illumination, simulations were performed with a pencil beam and a geometric focused beam [9,18].

Importance sampling and photon splitting methods
To increase the probability of backscattering events, the method of importance sampling was implemented and combined with photon splitting (see flowchart in Fig. 2) [24,34].Importance sampling is a probabilistic method to accelerate MC simulations for OCT by biasing random events leading to backscattering, thereby increasing the probability of photon detection and thus enhancing computational efficiency.An updated biased scattering function f B based on the HG probability density function f HG was then introduced.After a scattering event, this function modified the photon propagation direction to be oriented towards the detector, i.e., the collecting fiber system.The resulting bias was corrected by introducing a corresponding likelihood ratio as where L expresses the ratio between the probability of observing the biased scattering angle in an unbiased simulation θ S and in a biased simulation θ B .This likelihood ratio was tracked throughout all scattering events to ultimately adjust the weight of collected photon packets.Additional scattering, biased or unbiased, can be applied once the photon undergoes a first biased scattering event.To ensure the forward scattering of some photons, photon splitting was applied such that photons, who never experienced biased scattering, were split into two halves.One-half scattered in a biased manner towards the collecting fiber according to f B with a likelihood of L and the other unbiased by continuing the forward propagation pursuant to f HG with an assigned likelihood of L ′ = 1 − L. This last packet is considered to have never experienced biased scattering in Fig. 2 and would thus split again at the next scattering event.
Fig. 2. Flowchart describing the scattering process of the simulator.The randomly generated number is noted ξ i ∈ [0, 1] while p is the probability for biased versus unbiased scattering.

Optical configuration and tracking algorithm
Optical illumination and detection configurations were simulated using illumination and photon collection modules with identical properties of numerical aperture (NA), geometry, and spatial locations.Simulations of multiple A-lines were achieved by varying the lateral positions of the illumination and detection fibers.When the backscattered photon returned within the fiber's NA and radius, the photon was labeled as detected.The traveled distance of the photon packet in each different layer was tracked in addition to the fluence absorbed by the sample.This strategy allowed us to adjust the phase of the signal for media with different refractive indices afterward.Upon photon detection, the cumulative traveled distance (total pathlength TP), weight W, and likelihood L of each photon were saved.This method was implemented by modifying an existing open-source C-language software package that is available on the website of the Oregon Medical Laser Center [13].The MCOCT photon tracking algorithm flowchart is shown in Fig. 3 and describes the steps of a standard MC simulation (gray background) as well as our methodological contributions (white background).5 of the Appendix.

OCT signal reconstruction
For each A-line defined by a lateral (x) position of the illumination/detection fiber module, a number N of detected photons with their corresponding weight, likelihood, and total pathlength was acquired.These parameters were then used to reconstruct the OCT signal following a series of mathematical processes described in Fig. 4. To calculate the signal across the spectral bandwidth ∆λ at a central wavelength λ 0 , a discrete k array corresponding to the m th wavenumber (k m = 2πn/λ m ) was created with λ min = λ 0 − ∆λ 2 and λ max = λ 0 + ∆λ 2 .The collected signal S S (k m ) from the sample was then reconstructed as a function of wavelength and for all j detected photon packets using this expression: The reference arm signal amplitude S R (k m ) was calibrated as the average of the modulus | • | of the signal amplitude of the sample arm to ensure comparable power in both arms : By emulating a balanced detection configuration, the signal at the detector was generated by interfering the sample arm signals S S (k m ) with the reference arm signals S R (k m ) and therefore creating the interference fringe signal I(k m ) defined by Consequently, I(k m ) was defined as a function of the wavelength and required further Fourierdomain processing to generate the OCT signal as a function of depth by applying a Hanning window, zero padding, and discrete Fourier transform.These processes were performed in Matlab (Version R2021a, Natick, MA, USA).

Tissue-mimicking optical phantoms
Optical phantoms based on the optical properties of biological tissues were fabricated and characterized to validate this simulator.

Fabrication
To validate our simulator, optical phantoms mimicking biological tissue were developed with polydimethylsiloxane (PDMS) matrices.PDMS was selected for an accurate representation of the refractive index n ∼ 1.4 of biological tissue [39].The matrix demonstrates good durability and is visible with the OCT system at a wavelength of 1310 nm.Nigrosine was utilized as the absorbing agent and titanium dioxide powder as the scattering agent to mimic biological tissue (see Table 1).Three one-layer phantoms and one bi-layer phantom were fabricated with different optical properties by varying quantities of nigrosine and titanium dioxide.The bottom layer of the bi-layer phantom was created using the same recipe as the one-layer phantom #3, wherefore the properties of these layers are assumed equal.

Optical properties characterization
The absorption coefficient µ a and reduced scattering coefficient µ ′ s defined as µ ′ s = (1 − g)µ s were determined for the one-layer phantoms.To acquire the transmission and reflection spectra, a LAMBDA 1050 spectrophotometer from PerkinElmer (Waltham, MA, USA) was used.This commercial device is equipped with a single integration sphere capable of receiving both diffuse and specular lights.The collected spectra were processed using the inverse adding-doubling method to estimate µ a and µ ′ s [40,41].The determination of the anisotropy coefficient g was unattainable through direct measurement and was thus estimated by fitting OCT A-lines measured experimentally (see Sec. 2.4.3) to a theoretical model.In a first approximation, assuming a homogeneous turbid medium, the signal amplitude I(z) at a depth z can be modeled as a single exponential decay function following a variant of the Beer-Lambert law [42][43][44]: where I 0 is the incident light intensity.The factor 2 in the exponential decay accounts for the round-trip photon propagation typical to OCT.Assuming a direct measurement of µ a , the fitting procedure allowed to extract µ s .Visual inspection of the fit was performed in the range of depths z = [0.2,0.5] mm, where z = 0 referred to the sample's surface.This range was selected to avoid the impact of Fresnel's reflection.This model is also called the single scattering model, which assumes that the detected backscattered light undergoes only one scattering event in the sample.Measured optical properties are described in Table 2 and represent biological tissue for OCT imaging at 1310 nm [37,45].OCT imaging was performed with a VCSEL Fourier-domain system (Thorlabs, NJ, USA) at λ 0 = 1310 nm with ∆λ = 130 nm.Two lenses of different focal lengths were used: a short one with f = 30 mm and a long one with f = 75 mm.The position of the focus inside each optical phantom z f was determined separately for each focal lens and summarized in Table 3.For each experiment, OCT imaging generated 300 A-lines which were averaged for comparison with MCOCT.

MCOCT simulations for comparison with MCX simulations
To validate the beam focusing performance, the flux output (optical fluence rate) generated by MCOCT was compared to the one from MCX [18].The accumulated energy was stored within the 3D voxel matrix during the simulation.The 3D matrix with the deposited weight was normalized relative to the input power, resulting in the local absorption rate per watt (W) of the incident light in units of W/cm 3 /W incident .Dividing this rate by the absorption coefficients of the corresponding voxels resulted in the local fluence rate ϕ in W/cm 2 /W incident .The tissue was  represented in both simulators by a 200 × 200 × 300 voxel geometry with isotropic voxels of 0.001 cm 3 .The distance between the imaging lens and tissue surface z f ,I was set at 0.85 cm with a beam width at imaging lens D = 0.1 cm.In MCOCT, the illumination beam had a Gaussian profile at the imaging lens.In MCX, the trajectory was simulated as a hyperboloid of one sheet and the photon packets propagated in straight lines along the double-ruled surface, thereby mimicking a Gaussian-distributed photon launch [18].For both simulations, optical properties were µ a = 0.001 cm −1 , µ s = 0.01 cm −1 , n = 1, and g = 0.5.MCOCT and MCX simulations were executed using 10 6 photon packets and were run on a desktop computer equipped with an Intel Xeon E5 CPU at 3.50GHz, 16GB RAM and an NVIDIA TITAN Xp card.

MCOCT simulations for comparison with OCT imaging in phantoms
Optical properties from Table 2 were used to generate OCT A-lines with MCOCT.Importance sampling was incorporated with a bias coefficient of a = 0.9 and an additional bias probability of p = 0.8, as defined in [24,34].A minimum step s min = 10 −4 cm was used for these simulations.A scheme of 1000 equidistant illumination and detector pairs along the x-axis from x = −1 cm to x = 1 cm was used and resulted in simulating 1000 A-lines for a single B-scan.To minimize simulation time, a one-voxel slab of 2 × 2 × 0.5 cm 3 was used to model the one-layer phantoms.For comparison with the bi-layer phantom, the medium was modeled as a two-voxel slab of 2 × 2 × 1.8 cm 3 .For all simulations, the incident light beam diameter at the focusing lens D was set to 0.07 cm for imaging at λ 0 = 1310 nm.For OCT signal reconstruction, λ min = 1245 nm and λ max = 1375 nm.Optical properties were approximated to be constant over the OCT reconstruction spectrum.All simulations were executed using between 10 8 and 10 9 photon packets and were ran on the same desktop computer as above.Simulated A-lines are further averaged among detectors (1000) to reduce speckle contrast.

Results
Figure 5(a) and (b) show the optical fluence rates ϕ at y = 0 mm generated by MCOCT and MCX, respectively.In both simulations, the beam appears curved as expected for a Gaussian distribution.The simulation time was 60 seconds for MCOCT and one second using MCX. Figure 5 (c) shows a very good agreement between optical fluence rates simulated along the x-axis at the depth of the beam focus (z = 1.5 mm) by MCOCT (red) and MCX (blue).Figure 5 (d) presents the optical fluence rates simulated along the z-axis at a fixed x-position (x = 0 mm) by MCOCT (red) and MCX (blue).Both curves progress similarly around the beam waist but show more discrepancies with distance including a fastest decrease for MCX compared to MCOCT. Figure 6 and Fig. 7 show the comparison between A-lines generated by MCOCT (blue) and those measured in the 4 experimental phantoms (red) using an imaging lens with a focal length of f = 30 mm and f = 75 mm, respectively.At f = 30 mm, the attenuation of the MCOCT A-line follows a similar pattern to experimental data in the one-layer phantoms except in a region close to the surface of the sample where the signal intensity is lower in experimental data.In the bi-layer phantom (phantom #4, (d)), the change of optical properties from the superficial layer to the deep layer is marked by a transient change in signal intensity at a depth corresponding to the interface between layers (0.9 mm).A change in the attenuation slope is also observed with the change of optical properties between layers.For all comparisons, the accuracy of the simulated data also decreased with increasing depth as a result of lower signal intensity due to fewer collected photons.At f = 75 mm, signal intensity patterns were similar to that of f = 30 mm while signal intensity close to the surface presented a more accurate qualitative match.For both focal lengths, A-lines obtained by simulating a pencil beam and a geometric focused beam differed significantly from the A-line simulated with the Gaussian beam as well as from experimental OCT signals acquired in Phantom #1 (see Fig. 8 in Appendix).
Table 4 reports computation times for generating a B-scan of the phantoms presented in Fig. 6 and Fig. 7 along with the quantity of simulated and detected photons.Detected photons varied between 1.5 × 10 5 to 2.9 × 10 6 for a number of photons launched varying between 10 8 to 10 9 , yielding a ratio of detected/launched between 0.0525 and 0.33%.

Discussion
In this work, an OCT simulator (MCOCT) based on Monte Carlo methods was implemented and validated.The simulator was based on a Gaussian illumination profile and the importance sampling method as well as signal reconstruction in the Fourier domain.The numerical comparison with MCX [18] and experimental comparisons with OCT imaging in phantoms with known optical properties allowed us to assess and validate the performance of MCOCT.
In OCT systems, the high lateral resolution is achieved through tight focusing of the probing beam, which can be performed by the implementation of a Gaussian focusing profile.A method considering the wave nature of photons was developed by focusing the light along the curved trajectories in accordance with Gaussian optics to a diffraction-limited focal spot.In the simulation, a low scattering coefficient was selected to better visualize the curved Gaussian trajectory of photons.When compared to MCX hyperboloid [18], the diffraction-limited spot size at the focus w 0,n 2 matched well.Optical fluence rates from the two simulators showed an overlap in progression around the beam waist.With distance to the beam waist, the fluence rates varied, which may be due to the different approaches used for mimicking the Gaussian photon trajectory (hyperboloid vs. Gaussian optics).The fluence generated by MCOCT was similarly higher at the surface of the medium (z = 0 mm) and in-depth (z = 3 mm) compared to MCX.Using the importance sampling method in MCOCT reduced the number of photons that scattered in the forward direction, which also includes photons that should have exited the bottom boundary of the medium but were biased back to the surface.While this number is reduced for MCOCT, the fluence was corrected by the biased likelihood ratio and therefore does not explain the higher fluences observed with distance from the beam waist.If the likelihood ratios were not corrected for the biased scattering, the disparity between fluences calculated by MCOCT and MCX at the sample's surface would be different from that of the bottom.When focusing the light in MCOCT, the photon packet propagated by minimal steps allowing us to track the voxel in which the photon packet was located more accurately and adapt the direction in case of changing optical properties.This strategy allowed us to bypass the need to solve an equation of motion previously proposed [30,38].Then, only s min needs to be adapted and mainly depends on the voxel geometry.In MCOCT, the importance sampling method was used to increase the occurrence of backscattering, thereby reducing computational costs for OCT signal simulation, while in MCX all photons scattered according to the HG function.The difference in computation time between the two tools is mainly due to the computational architecture as MCX uses GPU resources while MCOCT uses CPU resources.The performance of the Gaussian illumination was also clear when comparing it with simulations performed with a pencil beam and a geometric focused beam.
When compared with OCT A-lines acquired in phantoms using a short focal length (30 mm), A-lines reconstructed by MCOCT differed within the first ∼0.2mm.This pattern was not observed at a longer focus (75 mm), where MCOCT matched OCT imaging in the four phantoms.This discrepancy is likely due to the reflection of light at the phantom surface, which is not considered in the simulator.When the surface has slight tilts, the coupling of reflected light with the incident light is higher for shorter focal lengths, leading to enhanced distortions in the A-line at the surface.It is also possible that a more accurate detection scheme [30] and scattering function [26,29] would provide a better match between MCOCT and experimental A-lines.The decrease in OCT signals with increasing penetration depth is expected as the probability of photon absorption also increases with greater propagation.A gradual deviation between A-lines generated by MCOCT and phantoms imaging occurred beyond a depth of about 1 to 1.5 mm.This difference is likely due to the fact that fewer photons reached these depths and were detected.The deviation was more prominent in one-layer phantoms imaged at the focal length of 30 mm compared to 75 mm.The reason behind this is unclear to us at present and is a current limitation of the simulator.For OCT applications in biological tissue, this limitation may not be a significant concern as typical OCT systems are designed to provide an imaging depth of approximately 1.5 mm [6].
The accurate characterization of signal intensity at interfaces is challenging when modeling complex geometry with multiple layers.When the change in refractive index occurs at the interface of two layers, photons undergo either reflection or refraction, which are mathematically described by Fresnel's equation [6].MCOCT signal intensity simulated proximal to the interface of the two layers in the bi-layer phantom did not match the one from OCT imaging at either focal lengths.Fresnel's reflection was not implemented in MCOCT and this discrepancy may be mitigated by considering such reflections as in previous works [30].
Minimizing computational costs is challenging when considering complex media with optical properties mimicking biological tissue.The computation time was higher for phantoms with higher scattering coefficients as these simulations presented more scattering events.For example, generating the OCT signal with the same quantity of detected photons for phantom #1 required higher computational resources than for phantom #3 as scattering was ∼6 times higher (see Table 4).Also, the quantity of photon packets launched when initiating the simulation dictates the computation time as well as the ratio of detected photons (0.0525 to 0.33%) that contributed to the OCT signal.As expected, simulations that used 10 9 photon packets (one-layer phantoms at f = 30 mm) were ∼10 times longer than simulations that launched 10 8 photon packets (one-layer phantoms at f = 75 mm).For the bi-layer phantom, a similar pattern was observed as the simulation using 4 × 10 8 photon packets (at f = 75 mm) took ∼4 times longer than one using 10 8 photon packets (at f = 30 mm).
In addition to the computational efficiency of the MC process, the number of voxels included in the geometry also affected the computational time of MCOCT.Each time a photon packet propagates by the distance prescribed by s min , MCOCT estimates its distance to the voxel boundary, which represents several occurrences when considering a geometry with a high number of voxels.This step is essential for geometry designed with multiple layers as optical properties have to be modified when propagating through different tissues.The voxel size also dictates the spatial resolution of the optical fluence rates.A strategy to minimize the number of crossing-boundary events is to use one voxel for each layer with the voxel size in the axial direction being equal to the layer thickness.This simple design can be used for simpler geometry such as rectangular slabs but is less appropriate for complex media.The axial resolution of the simulated OCT data is achieved through the imaging spectrum, whereas the quantity of simulated source/detector positions dictates the lateral resolution of the B-scan.The use of large voxels for simple geometry did not affect the lateral resolution of the OCT image.Another factor that influenced the computation time was the minimal step used to generate the Gaussian focusing beam.The selection of a short minimal step will increase the accuracy of individual photon trajectories but also increase the occurrences with which the trajectory has to be updated.Previous methods have been proposed to automatically adjust the step size during the simulation [38].In MCOCT, an adjustable step could be implemented but it is not clear how it would reduce computation time while maintaining high accuracy.Another time-consuming procedure in the Monte Carlo module was the random number generation that is needed for each scattering event.Generating one random number each time could be substituted with an a priori array of random numbers generated prior to the simulation.
Our work has several limitations.MCOCT simulations were run on a single CPU thread, which limited speed compared to GPU-based simulators [18,26,28,29].Accelerating methods such as parallel computing by Compute Unified Device Architecture (CUDA) will be explored for future versions of MCOCT.Using this architecture would also allow performing the initialization, simulation running, postprocessing, and visualization entirely in Matlab.Therefore, users would not need a separate UNIX-emulating console shell to run a compiler and the MC executable, nor would there be a need for intermediate input/output files, which additionally would save storage space.Parallelization would allow simulating a higher quantity of photon packets leading to higher accuracy of the OCT signals at higher sample depths.Noise and/or dispersion were not added to MCOCT signals before reconstruction, but could be added to the reconstructed OCT signal in future works.MCOCT used a scheme of 1000 illumination sources and detectors along the x-axis to generate a B-scan image.Since the geometry is described using a 3D voxel system, it is also possible to simulate a scan over the y-axis to reconstruct 3D volumes.However, a 3D scan would currently result in a significant increase in computation time.The demonstration of the ability to generate 3D volumes would be more relevant with a code based on parallelization using GPU.
Monte Carlo methods are the gold standard for studying photon propagation in turbid media [7,12,46] and have been instrumental in investigating the relationship between tissue optical properties and imaging results for several techniques including Fourier-domain OCT [22,26].An accurate analytical/numerical modeling of Fourier-domain OCT imaging is critical to study light-tissue interactions, including speckle contrast [47,48], dispersion [49,50], and multiple scattering [51,52].A better understanding of these optical phenomena via accurate simulations may help to design an instrument, optimize data acquisition specs and image reconstruction parameters.Numerical and in vitro experiments are often used to produce preliminary results that can guide the development of an in vivo experimental protocol to answer a scientific question.
While further enhancements remain necessary, our open-source simulator can be employed for model validation and system optimization by generating OCT signals across a wavelength spectrum to ascertain the appropriate NA to achieve the desired contrast and its observable depth.Simulating the impact of varying NA or different modes of a photonic lantern [53] can determine the feasibility of achieving a substantial contrast.Our tool can also contribute to the planning of in vivo experiments to provide preliminary estimates of specific biomarkers such as changes in retinal layer thickness to characterize particular physiological/pathological conditions [54].Previous applications [10,[19][20][21][22][23][24][25]27,28,30,31] and more recent ones involve the generation of OCT image datasets that can be used in the training of neural networks aimed at tackling the inverse problem in OCT.Some groups already proposed Monte Carlo simulation platforms to generate realistic OCT images for predictive model training [32,33].

Conclusion
In this work, a novel Monte Carlo simulator for Fourier-domain OCT (MCOCT) was introduced.The incorporation of a Gaussian light illumination scheme and the ability to increase backscattered event collection allowed us to achieve a more realistic approach to OCT imaging.Optical fluence rates generated by MCOCT were compared with an established MC simulator (MCX), which provided a numerical validation.OCT signals generated by MCOCT were also compared to OCT signals acquired at specific focal lengths in tissue-mimicking phantoms with known optical properties, which provided experimental validation.Future work involves leveraging GPU acceleration for enhanced computational efficiency and speed, as well as considering Fresnel's reflections proximal to interfaces in multi-layered geometries.MCOCT may further allow the designing of OCT system hardware and image processing for various medical and biomedical applications.

Appendix
Table 5. Variables used to describe Fig. 1 and Fig. 3 Table 5. Variables used to describe Fig. 1 and Fig.

Fig. 1 .
Fig. 1.Description of the geometry and parameters used to simulate the Gaussian focusing beam.The dashed lines show the trajectory in a crossing propagation.Variables are described in Table5of the Appendix.

Fig. 1 .
Fig. 1.Description of the geometry and parameters used to simulate the Gaussian focusing beam.The dashed lines show the trajectory in a crossing propagation.Variables are described in Table 5 of the Appendix.

Fig. 4 .
Fig. 4. Optical coherence tomography (OCT) signal reconstruction.Definitions of variables: TP N , total pathlength; W N , weight; L N , likelihood; N, N th collected photon packets; k, wavelength; S S (k m ), signal amplitude from the sample arm; S R (k m ), signal amplitude from the reference arm; I(k m ), reconstructed OCT signal; k m , m th wavenumber.

Fig. 5 .
Fig. 5. Optical fluence rates ϕ [W/cm 2 /W incident ] comparison between (a) MCOCT and (b) MCX simulated with the same optical properties.(c) Comparison of optical fluence rates along x at z = 1.5 mm and y = 0 mm generated by MCOCT (red) and MCX (blue).(d) Comparison of optical fluence rates along z at x = 0 mm and y = 0 mm simulated with MCOCT (red) and MCX (blue).

𝑛 1 ,Fig. 8 .
Fig. 8. Optical coherence tomography (OCT) A-line measured experimentally in Phantom #1 (red) compared to OCT A-lines simulated with a Gaussian (blue), geometric focused (orange), and pencil (green) beam using an imaging lens with a focal length (a) f = 30 and (b) f = 45 mm.