## Abstract

Optical coherence tomography (OCT) provides both structural and angiographic imaging modes. Because of its unique capabilities, OCT-based angiography has been increasingly adopted into small animal and human subject imaging. To support the development of the signal and image processing algorithms on which OCT-based angiography depends, we describe here a Monte Carlo-based model of the imaging approach. The model supports arbitrary three-dimensional vascular network geometries and incorporates methods to simulate OCT signal temporal decorrelation. With this model, it will be easier to compare the performance of existing and new angiographic signal processing algorithms, and to quantify the accuracy of vascular segmentation algorithms. The quantitative analysis of key algorithms within OCT-based angiography may, in turn, simplify the selection of algorithms in instrument design and accelerate the pace of new algorithm development.

© 2014 Optical Society of America

## 1. Introduction

The performance of optical coherence tomography (OCT)-based angiographic imaging has improved dramatically over the past 5 years [1, 2]. However, confounding artifacts and high noise levels continue to complicate the interpretation and segmentation of OCT angiographic datasets. Therefore, it is necessary to improve the signal and image processing algorithms upon which OCT-based angiography depends. To this end, both numerical and *in vivo* testbeds are useful. *In vivo* testbeds provide ultimate confirmation of performance, while numerical testbeds allow a more controlled evaluation of algorithms and optimization of their parameters. Numerical models based on Monte Carlo (MC) approaches have been used in the development of anatomical and polarization-sensitive OCT [3, 4]. However, these MC OCT models are inadequate for simulating microvascular signals for two reasons. First, OCT-based angiography is based on the detection of the time-varying nature of OCT signals, and the general approach of MC does not recapitulate these temporal variations. Second, OCT models are based on layered geometries that are too simplistic to model complex three-dimensional vascular networks.

In this paper, we describe a MC-based stimulation of OCT-based angiographic imaging of arbitrary three-dimensional vascular networks. The model supports multiple tissue types organized spatially into arbitrary three-dimensional geometries. We have also implemented Fourier-domain coherence gating directly into this model. As such, the model recapitulates many of the critical features of OCT-angiographic images including point-spread function broadening and vessel shadows. Most significantly, we have implemented MC-based simulation of OCT signal temporal variations within the regime wherein the flow-based displacement of scatterers is larger than the OCT system resolution (larger of transverse or axial). This allows modeling of angiographic signals. In the first two sections of this report, the methods for performing MC transport within arbitrary three-dimensional geometries (section 2) and the method for simulating both transverse and axial signal localization (section 3) are overviewed. In section 4, we describe the simulation of time-varying signals from the MC model. Finally, in section 5, we validate the approach through comparison of simulated images with datasets acquired *in vivo* in dorsal skinfold chamber preparations in mice.

## 2. Monte Carlo modeling of light transport within regions of arbitrary geometry

In prior works, MC models imposed a homogeneous or layered geometry [5, 6]. To support modeling of microvascular networks, it is necessary to extend these methods to support more complex region geometries. Our approach uses the same validated MC engine described in [5] but modifies the methods used to detect the intersection of a photon packet and a region boundary. First, we summarize the MC engine. Photon packets are launched from one surface location orthogonally into the sample. Photon packets begin with an initial weight, *W*, equal to 1. Once a photon packet is launched, the step length is selected randomly from an exponential distribution, the scale of which is defined by the optical properties of the tissue. The photon packet is then advanced along this direction orthogonal to the first surface until it either hits a boundary between regions, or it reaches the step length. If the photon packet reaches the step length, its weight is decreased by the ratio of the medium absorption to total interaction coefficient. The remaining photon packet is then scattered with a deflection angle, *θ*, and an azimuthal angle, *ψ*, which are statistically sampled from the Henyey-Greenstein scattering phase function and between 0 and 2*π*, respectively. If the photon packet intersects a region boundary along its advancement before reaching the step length, the packet is moved to the boundary. Here it is internally reflected or transmitted according a random number generator with probabilities matched to the Fresnel’s equation (including regimes for total internal reflection). After a transmission or an internal reflection, the photon packet completes the remaining dimensionless step distance divided by the new region absorption and scattering coefficients (if transmitted). A photon packet is discarded if it crosses an outside boundary or by a Russian Roulette random process once weight values fall below a predefined threshold limit. In the method of [5], regions are defined by z-axis locations and it is trivial to identify boundary transections.

To support a more complex geometry, we subdivide the sample into cuboidal voxels. We note that because these cuboidal voxels have boundaries that are always aligned to x/y axes, this method would not accurately estimate the specular (i.e., Fresnel) reflection from a tilted surface. In this approach, each voxel is assigned a specific tissue type with corresponding optical properties. The resolution of the geometry is defined by the voxel size, and can be reduced arbitrarily (at a computational cost) to simulate samples with higher resolution. For microvascular networks, the smallest vessels of which are on the order of 4–9 *μ*m, we adopt a voxel sizes of 3 *μ*m (27 *μ*m^{3}) [7, 8].

A boundary intersection algorithm was implemented to first identify the voxel boundaries that are intersected by an advancing photon packet. If there are no cube boundary intersections or if there are intersections but each intersection occurs between voxels associated with the same region (and therefore identical optical properties), the photon packet is transported the full step length. If there are cube boundary intersections within the photon packet path and one of these intersections occurs between two voxels of different regions, the photon packet is transported to the first intersection between differing voxels where it is internally reflected or transmitted. Analogously to the multilayered approach, after a transmission or an internal reflection, the photon packet continues to travel as described above the step distance adjusted by the new region total interaction coefficient.

To validate that the developed MC model provides similar results as validated layered codes [5], we compared the absorption distribution of light in a three-dimensional medium. To record absorption in the cross-sectional plane *x* = 0, a *k* × *l* grid was defined where each element has a dimension of 5 × 5 *μ*m. Photon packets were launched from one surface location at the origin of the medium and absorption was scored during propagation according to the appropriate grid element in a two-dimensional matrix *A*(*y*, *z*), where 1 ≤ *y* ≤ *k* and 1 ≤ *z* ≤ *l*. Figures 1(a) and 1(c) show the geometry of the models used to simulate the absorption distribution, and Figs. 1(b) and 1(d) show the absorption distribution for the multilayered MC code [5] and the newly developed code, respectively. Note that cylindrical regions were included within the new model to exercise the boundary condition code for arbitrary geometries, but that the optical properties of each cylindrical region are identical to its surrounding medium and as such the resulting absorption distribution should match that calculated by the multilayered MC code.

## 3. Modeling OCT transverse and axial signal localization

The aforementioned methods allow MC-based calculation of light transport in arbitrary three-dimensional structures. To extend this model to simulate OCT imaging, it is necessary to include methods for proving both transverse and axial signal localization (e.g., the imaging component of OCT). Transverse resolution can be incorporated as has been done in prior models [3, 9]. Briefly, photon packets are launched with a location defined by a uniform circular distribution. The angular position of photon packets is sampled uniformly between 0 and 2*π* while the radial position is sampled between 0 and a 10-*μ*m beam radius with a square root distribution. On collection, exiting photon packets are included only if they meet the spatial and angular constraints of a detection aperture. In our model, a spatial filter with a 10-*μ*m radius (uniform intensity) and an angular filter limited to a +/− 5° detection angle [3] was adopted. The location of the launch and collection apertures is translated across the imaging surface to simulate beam scanning.

In our approach, we have modeled directly the generation of wavelength-dependent fringe signals, and we then use Fourier analysis to convert these fringes to depth-resolved A-lines. The simulation tool records the square root of the diffuse reflectance (weight) and the total path length of each photon packet that passes the detection aperture. At the conclusion of the MC simulation for a given beam location (e.g, the modeling of one A-line), we obtain a set of detected photons packets, each with a diffuse reflectance amplitude and a path length. We then construct the received sample signal (amplitude and phase) as a function of wavelength according to

where*R*and

_{n}*l*respectively represent the recorded diffuse reflectance and path length of the

_{n}*n*

^{th}detected photon packet.

*k*corresponds to the wavenumber and

*G*(

*k*), to the Gaussian light source spectrum defined as where

*k*

_{0}represents the center wavenumber of the light source and

*dk*the spectral bandwidth. The same diffuse reflectance,

*R*, and path length,

_{n}*l*, arrays are used for each

_{n}*k*. We calculate

*S*across a discrete

*k*array between

*k*

_{start}(

*λ*= 1220 nm) and

*k*

_{stop}(

*λ*= 1380 nm) with equal spacing in

*k*.

The reference field is calculated as

where*α*is a parameter chosen to scale the total sample arm power (|

*S*(

*k*)

^{2}|) between 0.001–0.01% of the reference arm power (|

*R*(

*k*)|

^{2}). The value of the parameter

*α*is constant across a simulation for all A-lines.

We interfere the sample and reference arm fields in a balanced detection configuration to generate a fringe signal given by

We then follow conventional Fourier-domain processing by applying a window (Hamming) and Fourier transform to give the OCT signal as a function of depth. Figure 2(a) shows an OCT interference signal, *I*_{D}(*k*), as a function of wavelength obtained with a two-layered medium (air and cutaneous tissue), where the first layer starts at 0 *μ*m and the second, at 50 *μ*m. Figure 2(b) presents the corresponding *i*_{D}(*z*) signal as a function of depth following Fourier transform. Inspection of Fig. 2(b) reveals a typical OCT signal where the two regions are clearly visible.

To validate the OCT signal localization methods, we generated simulated OCT structural images of a simple phantom structure with a homogeneous background and two vessel inclusions. Figure 3(a) shows the geometry of the medium comprising the two blood vessels with a diameter of 150 *μ*m and 100 *μ*m. Figure 3(b) presents the corresponding structural OCT image showing the vessels with the typical shadowing effect and omnipresent speckle noise.

## 4. Simulating blood flow

OCT angiography is based on the acquisition of repeated measurements from the same location and the analysis of these measurements to identify time-varying signals that result from flowing blood. To simulate angiographic OCT images, it is therefore necessary to derive time-series measurements that are constant or modulated according to their degree of interaction with intravascular scatterers. In this work, we make an assumption that the intravascular scatterers have translated by more than the larger of the axial or transverse resolutions. As a result there is no correlation between the OCT signals scattered from vessels across two time points. This accurately describes many of the angiographic OCT approaches that employ large time separations between measurements to achieve sensitivity for slower flowing vessels.

To model temporal decorrelation of intravascular scatterers, we used the following approach. First, during MC light transport, we flag each photon packet (using the vessel flag variable, *F _{n}*) that scatters within a voxel associated with the region defining the intravascular space (section 2). Then, during calculation of the OCT signal

*I*

_{D1}(

*k*), we include in each photon packet that is flagged a random phasor,

*ϕ*

_{1n}, between 0 and 2

*π*within the expontential. To simulate a second measurement,

*I*

_{D2}(

*k*), at the same location, the same MC solution (

*R*and

_{n}*l*) is used but a new set of random phasors,

_{n}*ϕ*

_{2n}, are applied. In this way, the OCT signals derived from photon packets that have not interacted with intravascular scatterers are constant in time, but those signals that are derived, at least in part, from photon packets scattered from intravascular scatterers are proportionally modulated. In this work, the OCT angiography contrast is then calculated from the depth-directional variance of the differential signal between

*i*

_{D1}(

*z*) and

*i*

_{D2}(

*z*) [2]. However, it is also possible to model other angiographic algorithms that operate within the regime of full signal decorrelation between measurements. We note that the effect of vessel size on the OCT signal is accounted for by the inherent statistical nature of Monte Carlo methods and not, for example, by adjusting the magnitude of the phase (

*ϕ*) in the random phasor according to vessel size. A smaller vessel will have a larger proportion of photon packets that traversed the vessel without a scattering (and are therefore not phase modulated) than would a larger vessel. We also note that the simulation of time-series data does not require that the MC simulation of light transport be recalculated, and thus adds minimal additional computational burden.

_{n}The MC model algorithm is summarized in Fig. 4. For clarity, the function of each step is represented by its box line style. Boxes outlined by a full line describe photon transport in current MC models [5], by a dashed line, photon propagation in a discretized medium, by a dotted line, OCT-based signal localization, and by a dash-dotted line, modeling of flowing blood. This algorithm was implemented in a C language environment and was developed to provide acceleration of simulation time by offering the option to run on multithreaded Central Processing Units (CPU) or CUDA Graphics Processing Units (GPU) [10] depending on available hardware. The use of parallel threads enables many photon packets to be simulated simultaneously. For instance, simulation time to generate Fig. 1(b) on an Intel Core i7 processor with the existing downloadable model was 152 seconds with a single thread, while simulation time to obtain Fig. 1(d) with the new model required 3 seconds with 12 threads. The improvements in algorithm speed are necessary to make complex three dimensional simulations practical.

While simulation time is largely dependent upon hardware, it is also affected by the number of launched photon packets and cuboids implemented in the model. While the use of a large number of photons and cuboids in simulations improves signal precision and three-dimensional shape resolution, it may rapidly become computationally impractical. The user must therefore appropriately choose optimal parameters for a specific application to obtain desired data within a reasonable timeframe.

## 5. Validation

To validate the model, we compared the appearence and signal statistics of a simulated angiographic image with that of angiographic images acquired *in vivo* (Figs. 5(d) and 5(h)). The simulated phantom included vessels designed to match that of a reference *in vivo* image acquired with an existing OCT-based system on a mouse skin tumor through a dorsal skinfold window preparation. The *in vivo* data was acquired at an A-line rate of 50 kHz with a 1300 nm system. A smaller region of the OCT image containing two larger blood vessels (diameters of approximately 250 *μ*m (left) and 70 *μ*m) was extracted.

Figure 5 shows cross-sectional images of the two blood vessels simulated using different optical properties based on published values [12, 13]. While values for absorption and scattering coefficients are in good agreement for wavelengths above ∼ 1,050 nm [13], the anisotropy factor in [12] remains larger (*g* ≃ 0.995 [12] compared to *g* ≃ 0.97 [13]). Figures 5(a) and 5(e) were generated based on values published in [13] while Figs. 5(c) and 5(g) were obtained from data presented in [12]. Figures 5(b) and 5(f) were simulated using an intermediate anisotropy value. All results show the two blood vessels with the OCT shadowing effect below each vessel. Based on these simulations, the scattering coefficient for blood, *μ*_{s}, was set to 650 cm^{−1}, the absorption coefficient, *μ*_{a}, to 5 cm^{−1}, and the anisotropy factor, *g*, to 0.9888.

Figure 6 compares the angiographic signal statistics in simulated and *in vivo* datasets. Histograms were generated over a region of interest (ROI) in Figs. 6(a) and 6(b) with an axial depth of 100 *μ*m and a width of 200 *μ*m. Inspection reveals agreement in signal distributions (Figs. 6(c) and 6(d)).

## 6. Discussion and conclusion

This paper describes a new MC-based simulation tool for computing OCT angiography images. The tool offers three key innovations relative to existing numerical models of OCT. First, a voxel-based region definition allows arbitrary three-dimensional structures including complex vascular networks. Second, MC methods are used to capture photon packets and their propagation distance, and Fourier-domain coherence gating methods are used to generate depth-resolved fringes. Third, a method allowing simulation of flow-induced temporal decorrelation was implemented which enables angiographic image generation.

While this tool provides simulation of angiographic signals, it does not accurately recapitulate the noise present in angiographic imaging of static scatterers. This will be a focus of future work, and would enable direct comparison of the contrast of angiographic signals relative to background for varying angiographic algorithms. The tool could be further optimized to include the Gaussian beam profile on the launched and collection apertures; the current system uses a circular aperture. While these changes to the aperture will allow a more accurate modeling of the system transverse resolution, they will not affect the statistics of the angiographic signals. This tool can be complemented by adding methods to mathematically generate realistic microvascular networks for simulation of three-dimensional OCT imaging. This would allow segmentation algorithms to be optimized on datasets in which the true vessel geometries are known *a priori*. Finally, this model may be extendable to allow the modeling of partial temporal decorrelation, and would allow the testing and development of quantitative flow measurement algorithms in addition to the morphological angiography described in this work.

## Acknowledgments

This project was supported by the Center for Biomedical OCT Research and Translation through Grant Number P41EB015903, awarded by the National Center for Research Resources and the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health. This work was also supported by NIH grant R01CA163528 and by the Natural Sciences and Engineering Research Council of Canada.

## References and links

**1. **B. J. Vakoc, D. Fukumura, R. K. Jain, and B. E. Bouma, “Cancer imaging by optical coherence tomography: preclinical progress and clinical potential,” Nat. Rev. Cancer **12**(5), 363–368 (2012). [CrossRef] [PubMed]

**2. **B. J. Vakoc, R. M. Lanning, J. A. Tyrrell, T. P. Padera, L. A. Bartlett, T. Stylianopoulos, L. L. Munn, G. J. Tearney, D. Fukumura, R. K. Jain, and B. E. Bouma, “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nat. Med. **15**(10), 1219–1223 (2009). [CrossRef] [PubMed]

**3. **G. Yao and L. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol. **44**(9), 2307–2320 (1999). [CrossRef] [PubMed]

**4. **I. Meglinski, M. Kirillin, V. Kuzmin, and R. Myllyla, “Simulation of polarization-sensitive optical coherence tomography images by a Monte Carlo method,” Opt. Lett. **33**(14), 1581–1583 (2008) [CrossRef] [PubMed]

**5. **L. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. **47**(2), 131–146 (1995). [CrossRef] [PubMed]

**6. **L. Wang and H.-I. Wu, *Biomedical Optics: Principles and Imaging* (Wiley, 2007).

**7. **R. A. Freitas Jr., *Nanomedicine, Vol. I: Basic Capabilities* (Landes Bioscience, 1999).

**8. **M. P. Wiedeman, “Lengths and diameters of peripheral arterial vessels in the living animal,” Circ. Reach. **10**, 686–690 (1962). [CrossRef]

**9. **M. Y. Kirillin, G. Parhat, E. A. Sergeeva, M. C. Kolios, and A. Vitkin, “Speckle statistics in OCT images: Monte Carlo simulations and experimental studies,” Opt. Lett. **39**(12), 3472–3475 (2014) [CrossRef] [PubMed]

**10. **A. Doronin and I. Meglinski, “Online object oriented Monte Carlo computational tool for the needs of biomedical optics,” Biomed. Opt. Express **2**(9), 2461–2469 (2011). [CrossRef] [PubMed]

**11. **A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: a review,” J. Innov. Opt. Health Sci. **4**(1), 9–38 (2011). [CrossRef]

**12. **A. N. Yaroslavsky, I. V. Yaroslavsky, T. Goldbach, and H.-J. Schwarzmair, “The optical properties of blood in the near infrared spectral range,” Proc. SPIE **2678**, 314–332 (1996). [CrossRef]

**13. **M. Friebel, A. Roggan, G. Muller, and M. Meinke, “Determination of optical properties of human blood in the spectral range of 250 to 1100 nm using Monte Carlo simulations with hematocrit-dependent effective scattering phase functions,” J. Biomed. Opt. **11**(3), 034021 (2006). [CrossRef]