Monte Carlo simulation studies of optical coherence tomography (OCT) and optical Doppler tomography (ODT)

A Monte Carlo simulation model has been developed to study how OCT and ODT imaging results are influenced by multiple scattering effects. For shallow optical depths [less than 3 mean free path (mfp) units] in intralipid (g equals 0.7), the number of detected backscattered photons was found to follow the extinction-single-backscatter model, and accurate mean values and small standard deviations for axial and lateral point spread functions were observed. For deeper depths, the positional accuracy and precision of backscatter was rapidly lost, and the number of detected photons no longer decreased with increasing focus depth in the nonabsorbing medium. For strongly forward-directed scattering in blood (g equals 0.99), depth profiles of the number of detected backscattered photons only approached the extinction-single- backscatter model for very shallow depths (optical depth less than 2 mfp units, metric depth less than 13 micrometers). However, quite accurate and precise point spread functions were observed even for large optical depths (40 mfp units). Doppler detection of blood flow was studied by simulating a 100 micrometer diameter horizontal blood vessel placed 200 micrometers below the surface in 2% intralipid. Simulated depth profiles of average Doppler frequency demonstrated good accuracy in absolute velocity values and localization of flow boundaries (3 - 4% deviation). Doppler frequency noise was observed in backscattering from regions underneath the vessel.


INTRODUCTION
Signal attenuation and loss of resolution due to multiple scattering severely limit the depth penetration of OCT, and make it difficult to image deep into strongly scattering biological tissues.Mathematical models have been developed in an attempt to increase the understanding of the OCT imaging process and thereby possibly enable development of better imaging instrumentation and data processing procedures.The extinction-single-backscatter model represents the OCT imaging process as an incoming beam which is continuously attenuated according to the total extinction coefficient of the propagation medium'5.After a single backscattering event, light is backpropagated with continuous attenuation towards the detector.Attempts have been made to include multiple scattering effects into such analytical models6, but meaningful comparison with experiment remains difficult since only simple geometries can be considered.For modeling of more realistic geometries, Monte Carlo simulation represents the only feasible alternative.
Monte Carlo modeling is also a suitable approach for studying the effects of multiple scattering on Doppler flow measurements implemented in low coherence reflectometry.This experimental technique is known as optical Doppler tomography and has yielded surprisingly good results7'2, considering that multiple scattering effects are bound to influence the measurement even from the smallest blood vessels, since the scattering mean free path in blood is only about 7 .tm.It is not difficult to construct examples from nonuniform velocity fields such as the parabolic velocity field for laminar flow in cylindrical vessels, where the total Doppler shift reported by a photon that undergoes two scattering events, will be very different from the exact value reported from a single backscattering event.Furthermore, the apparent backscatter location which will be ascribed to the detected photon, will be different from either ofthe two scattering locations.
In the present report we show some of the results obtained by a Monte Carlo model developed to study both OCT and ODT'3' ' Full 3D stochastic simulation was implemented in the medium.The OCT detector was modeled with low coherence interferometric path length gating and confocal selectivity determined by the antenna theorem for a heterodyne receiver, which limits the effective size of the OCT detector to twice the diameter of an ideal confocal detectors ' 15 This represents nearly full confocality, and OCT has been said to perform as a confocal microscope with an additional narrow time gate -'5O fs) to further increase the selectivity ofphoton detection5.Output data from the Monte Carlo simulation were the number of detected photons, as well as their accumulated path lengths and Doppler shifts.Also included were complete records of spatial locations for all scattering interactions which allow calculation of point spread functions for backscattering localization relative to the nominal focus position of the probe.With this model we can study how instrumentation parameters and properties of the flowing medium and its surroundings might influence the positional accuracy of OCT and ODT measurements and the accuracy of Doppler velocity determinations.

Probe geometry
A geometric model of the combined emission and detection probe of a fiberoptic OCT/ODT instrument is shown in Figure 1.The probe is pointed at an incidence angle of a (typically 15°) since we want to 2r simulate Doppler measurements of blood flow in a horizontal tube immersed at an axial depth of zo in a scattering medium.The probe is aimed in the direction of point P on the z-axis, but due to refraction at the surface of the medium, the focus point is shifted to point ' at depth z' where z'=kz.The constant k is determined by using Snell's law for the refraction of light propagating into a medium of relative index of refraction n (=1 .37).The probe has a numerical aperture NA=sinO, defined by the divergence of the light from the fiber core F and the focusing lens L. The surfaces Sj and S2 defme the spherical wavefronts of the beam diverging from the fiber tip and being focused by the lens L to converge towards point P (in air).Simulated photons are launched from the probe in such a way that the numerical aperture is filled with a Gaussian intensity profile.The probe acts as a confocal detector, since the antenna theorem for a heterodyne receiver limits the effective size of an interferometric receiver to twice the diameter of a confocal detectors' 15 Figure 1 .Geometry of probe and medium.This feature is implemented through the confocality angle 8. Photons that are backscattered from the medium therefore have to meet three criteria to be detected: (1) intersect the surface S2 defmed by the numerical aperture and distance from point P', (2) deviate from normal incidence on S2 by less than the confocality angle 5, and (3) have a total optical path length that deviates by less than the source coherence length L, from the nominal pathlength for a photon launched from S2 and scattered from P' back to S2 in a single (back)scattering event.Intralipid was chosen as a model medium that could easily be used experimentally to verify predictions resulting from the Monte Carlo study.Both intralipid and blood were modeled as homogeneous media, each characterised by three optical parameters, i.e. the absorption coefficient the scattering coefficient p, and the anisotropy parameter g.Based on limited available literature data for 2=850 nm, the following values were chosen'3' ' for I % intralipid: p=O, p=2.O mm1, g=O.7; and for blood: p,=O.75 mm' , fls=150 mm', g=O.99.

Generating photon histories
Based on a chosen orientation and numerical aperture of the input probe (usually NA=O.2), a set of photon histories was simulated for each focus position of the probe.Each photon launched from the probe was propagated into and within the medium while keeping account of its accumulating pathlength in medium-equivalent physical distance units.Within the medium, the length of each path segment was determined by stochastically sampling the path-length distribution according to the method of Buslenko et a116 as described previously'7.At each scattering interaction with the medium, the scattering angle was determined by stochastically sampling the Henyey-Greenstein phase function'8' 19 For each scattering interaction, the position of the scattering event was recorded, and the accumulated path-length was incremented by the length of the latest path segment.At each scattering interaction with the moving medium in the vessel the resulting Doppler shift was determined as: (1) where /; and are the incoming and scattered wave vectors for this particular interaction, respectively, and i is the velocity ofthe medium at the position ofthe interaction.
For the generation of photon histories, maximum values for detection parameters such as numerical aperture (NA), source coherence length (La) 1fld confocality angle ( were pre-set by the user.This made it possible to define for each depth position a minimum and maximum path length for backscattered photons that may be detected.Propagation of any photon was terminated and its history discarded as soon as it could not make it back to the surface S2 without exceeding the maximum path-length for the depth in question. Backscattered photons that emerged through the surface of the medium were refracted and propagated to see if they hit the surface S2 defmed by the maximum numerical aperture of the detector.The incidence angle of the photon at S2 was tested to see if it was within the maximum confocality angle 8m• Path histories were written to disk file for each photon that satisfied the set of maximum values set for the detection parameters NA, L and 8.Each history included (x,y,z' )-positions and possible Doppler shift for all scattering events, spatial and angular coordinates of the photon's intersection with the detector surface S2 , and the accumulated Doppler shift as well as the accumulated path length.The accumulated Doppler shift from all flj scattering interactions in the medium was defmed as the Doppler frequency, COD, ofthe detected photon:

Analyzing photon histories
The data processing program was implemented with the AVS image analysis software (Advanced Visualisation System, Waltham, MA) in the form of a number of specially-developed AVS modules that integrate into the total AVS software package and utilise standard AVS modules for data processing and graphic presentation.Special-purpose modules were written to: read photon histories from disk file for a specified subset, or all, of the focus positions simulated in a specified run of the photon histories generator.
determine which of the read photons satisfy the present detection criteria, based on interactively specified values for the detection parameters NA, L and 6 count the number of photons detected for each (x,z') position in a scan series.
generate a histogram ofthe number of scattering events for photons detected for the specified focus position (x,zD. -display graphically the ensemble of individual path segments for the detected photons from one (x,z' ) focus position of the probe.This module utilises other standard AVS graphics modules that allow scaling and rotation of the 3D output graph.
display for all detected photons from a certain focus position (x,z'), all the interaction points that the photons have had in the medium.Alternatively, the module displays only the interaction point at maximum depth for each individual photon.These are the points of backscattering, and this output is defined as the point spread function.The geometry can be viewed as an xz'-image with the data collapsed along the y-axis, or as a yz'-image with the data collapsed along the x-axis.For each of these two images, a collapsed one-dimensional distribution corresponding to the axial (zPSF ) or lateral point spread function (xPSF or yPSF, respectively) can be generated.
generate the Doppler frequency distribution, which for a specified position (x,z' ), shows a histogram of the Doppler frequencies co, for the photons detected at that position.
calculate the average Doppler frequency for photons detected at various (x,z' )-positions in a z'-scan.The average is computed as: WD0'DiDj (3) ni =i where flj i5 the number of photons detected for each position (x,z').

Signal attenuation and localization
Attenuation and localization of signal was studied by aiming the probe into the homogeneous intralipid or blood medium, i.e. with the blood vessel shown in Figure 1 removed.The probe was aimed at an incidence angle of 15° in order to correspond to situations where Doppler shifts also would be determined.
Figure 2 shows intensity profiles as a function of focus depth in intralipid and blood.Also shown are curves determined by the extinction-single-backscatter model, i.e. according to the equation: exp(-2ptz '/coslS 9(broken lines) wherep1=p+p is the total attenuation coefficient.For intralipid the simulated data are well represented by the extinction-single-backscatter model down to an optical depth of about 3 mean free path (mfp) units (i.e.ptZ ' 3), but for greater depths, the number of detected photons levels offand shows essentially no further decrease with increasing depth in the nonabsorbing medium.
The transition between the region where the data are well represented by the extinction model and the region of constant backscatter intensity for greater depths, was found to depend strongly on the confocality angle 8.As S was decreased, the data followed the extinction model to increasingly larger depths (results not shown).A smaller value of S represents a higher selectivity in photon detection, in that photons in order to be detected must have a direction that deviates less than the limiting value from the ideal normal to the spherical wavefronts S2 and Si at the lens of the probe (Fig. 1).Setting a smaller confocality angle 5 therefore causes only minimally scattered photons to be detected.Such photons have followed tracks that deviate only minimally from the track of a photon experiencing only a single (back)scattering event at the nominal focus depth.Thus, minimally scattered photons have suffered only small changes of direction and/or have only short path segments that deviate significantly from the straight track represented by a single (back)scattering event at the nominal focus depth.

Focus depth (mfp)
Figure 2. Relative number of photons detected from intralipid and blood as a function of optical depth (p8Z') in mfp units.
The increasingly more selective detection of only minimally scattered photons that was achieved by decreasing the confocality angle 5, caused a corresponding decrease in the number of diffuse photons detected from depths beyond the valid range of the extinction model.Due to run time limitations, it was not possible to simulate a detector having the full selectivity of an ideal heterodyne receiver, which would correspond to confocality of 28=0.05°foraNA=O.2 detector with a 3 mm working distance.
Other detection parameters such as the numerical aperture (NA) of the detector and the coherence length of the light source (La) also influenced the number of detected photons, but only in a way that shifted curves like those in Figure 2 along the ordinate axis, i.e.only changed the number of detected photons by the same factor for all optical depths.Data such as those in Figure 2 were therefore generated using the maximum values of NA=O.4 and L=5O im in order to have the highest possible number of detected photons.
The simulated data for blood show a decrease in detected photons with increasing depth that is not at all well represented by the extinction-single-backscatter model, not even for very shallow depths.Decreasing the confocality angle 8 also in this case increased the initial steepness of the curve (results not shown) but not to the extent that the initial slope came close to that of the extinction model.Thus, for the very strongly forward-directed scattering in blood (g=O.99), the selectivity of detection could not be increased by decreasing S to the extent that only minimally scattered photons were detected.
Consequently, it was not possible for blood to distinguish between detection of minimally scattered photons at low optical depths and detection ofdiffusely scattered photons at greater optical depths.
Figure 3 shows individual tracks of detected backscattered photons from intralipid and blood for various optical depths.In intralipid, panel A shows well-localised backscattering from the focus position at an optical depth of 3 mfp units (300 im in 5 % intralipid, corresponding to reduced optical depth of 0.9), even if a few of the detected photons reveal tracks that deviate considerably from the path of single-(back)scattering events.As the depth in intralipid is increased to 7 mfp units (panel C, 700 .tm in 5 % intralipid, corresponding to reduced optical depth of 2.1), all detected photons have very jagged The highly forward-directed scattering in blood (g0.99) preserves directionality of the photon beam through multiple scattering events.Panel B represents an optical depth of 22.5 (150 tm depth in blood, corresponding to reduced optical depth of 0.2) and shows localization of backscattering that is nearly as good as for an optical depth of 3 in intralipid (panel A).At an optical depth of 105 (panel D), corresponding to a reduced optical depth of I .0,localization of backscatter signal is about to be completely lost.
For data such as those shown in Figure 3, point spread functions (PSFs) were generated that show the distribution of backscattering positions for detected photons in either of the three coordinate directions.The graphs in Figure 4 represent mean values of zPSF and yPSF as a function of optical depth in intralipid and blood.Also shown is the expected mean value of zPSF (broken line) if backscattering always took place at the nominal focus point P' (Fig. 1).

A C D
In intralipid, localization of backscatter is only preserved down to an optical depth of about 3 mfp units, corresponding to the detection of only minimally scattered photons, as explained for Figure 2. At deeper focus depths in intralipid photons are backscattered from considerably more shallow depths than the nominal focus position.This is evidently because multiple scattering results in very jagged tracks in intralipid.Since only photons with the nominally correct accumulated path length are accepted, detected photons have not penetrated to the full focus depth.Error bars in Figure 4 represent 1 standard deviation of the distributions of backscatter locations.The relative standard deviation (i.e.precision) ofthe zPSF increases rapidly when diffusely scattered photons start to be detected at optical depths exceeding 3 mfp units.For optical depths of 20 mfp units, the standard deviation of detected backscatter position has increased to about 30% of the mean depth of actual backscatter positions.The mean backscatter position along the ydirection (see Fig. 1) is zero as long as only minimally scattered photons are detected, but assumes a gradually more positive value for greater depths.This is due to the oblique angle of incidence of the probe.Thus, photons interacting with the medium in the region corresponding to positive y-coordinates, have a larger probability of being detected than photons scattered into the half-space corresponding to negative y-coordinates.The standard deviations of the yPSF are about 50 % larger than the corresponding standard deviation for zPSF.A reasonable explanation for this is that the gating on total pathlength to be within of the nominal path for the chosen focus depth, limits the range of z-positions of backscattering events more efficiently than it limits the spread of lateral positions.
The mean of the xPSF remained at zero for all focus depths, as expected for symmetry reasons, and the standard deviations ofxPSF and yPSF were comparable (results not shown).
In contrast to intralipid, backscattering positions in blood were well localised to the nominal focus position even for optical depths greater than 35 mfp units.The detected backscattering depth was only about 3 % shorter than the nominal depth and had standard deviations of about 3 % at a depth of 35 mfp units.Figure 5 shows average Doppler frequencies according to eq. 3 for photons detected as the focus of the probe was moved along the z-axis from a nominal depth of 190 .tmthrough 390 .tm.The velocity distribution in the vessel was taken to be the parabolic profile for laminar flow in a cylindrical tube, with a maximum velocity of 2 mm/s.The broken curve in each panel shows the Doppler frequency distribution expected from the specified parabolic velocity profile.According to eq. 1 the detected Doppler frequency for a backscattering event where Ak 1E, -k -2E1 is proportional to cos(90 °-a), and therefore increases as the incidence angle of the probe is increased.The average Doppler frequency profiles shown in Figure 5 closely followed the expected parabolic profile and demonstrated the expected increase with increasing incidence angle ofthe probe.
Photons backscattered from (static) intralipid underneath the vessel showed random noise Doppler frequencies.Both positive and negative values were found, and the noise level seemed to be about the same for the two different incidence angles.
Figure 6 shows standard deviations of Doppler frequency distributions determined for photons detected from various depths through and under the vessel.Values are shown for two different numerical apertures of the detector.The noise level from underneath the vessel is of the order of 400 to 600 Hz and seems to be essentially independent of incidence angle.The noise values furthermore seem to be the same for the two different numerical apertures of the detector, even if values for the lower numerical aperture show larger random fluctuations due to fewer photons detected.Standard deviations determined from the flow region seemed to follow the parabolic velocity profile and were significantly larger for larger numerical apertures of the detector (circles versus triangles, Fig. 6).The latter finding can be explained by reference to Figure 1 , where it will be seen that the Doppler shift resulting from a backscattering event at the focus position will vary between 2kv .cos(90°-a-O) and 2kv cos(90 °-a+ 0), i.e. show larger variation for the larger 0 corresponding to higher numerical aperture.Interestingly, the standard deviation of Doppler noise from underneath the vessel is larger than the standard deviation of Doppler frequencies detected from the flow region.

DISCUSSION
The present study was based on modeling of homogeneous media characterised only by absorption and scattering coefficients in addition to the anisotropy parameter which defmes the degree of forward-directed scattering.Monte Carlo simulation of photon propagation and Doppler effects within the medium were implemented in a similar way as used for studies of laser Doppler flowmetry2022.This was combined with a geometrical optics model of the probe and interferometric detection which takes into consideration the antenna theorem of a heterodyne receivers' 15The interferometric gating on path length and the high degree of confocality inherent in the heterodyne receiver constitute a highly selective detector in terms of accepting only photons backscattered from the focus region of the probe.
Our simulation results (Fig. 2) show that for a moderately anisotropic medium such as intralipid (g=O.7), the highly selective confocal-like detection inherent in OCT measurements causes only minimally scattered photons to be detected from depths down to several scattering mean free paths.Within this region of minimal scattering, the detected intensity follows the simple extinction model exp(-2p,z'/cosa).Minimally scattered photons have fewer than average number of scattering events for the optical depth in question, have more forward-directed scattering than average, and have only short The highly forward-directed scattering in blood combined with the difficulty in studying very high confocality due to low photon yield, make it difficult to identify a region of minimal scattering in blood where detected intensity is governed by the extinction model.However, the highly forward-directed scattering in blood preserves localization of backscatter signal to the nominal focus position even for great optical depths.Since directionality is preserved in spite of multiple scattering, the depth resolution may be improved by a short source coherence length, and accurate position localization is obtained down to considerable optical depth (Fig. 3, 4).
Loss of directionality in photon propagation is more related to the reduced optical depth pz '(l-g) than to the optical depth per se.An arbitrarily chosen criterion of a 10 % relative standard deviation (CV) in the depth point spread function was found to correspond to reduced optical depths of 0.9 and 0.6 units of reduced scattering mean free paths in intralipid and blood, respectively (results not shown).This indicates that loss of depth resolution essentially occurs at a reduced optical depth of 1 mfp unit.The optical depths corresponding to 1 unit reduced optical depth are, however, widely different for intralipid and blood, i.e. 3 .3for intralipid corresponding to 0.8 mm depth in 2 % intralipid, and 100 in blood corresponding to 0.7 mm depth.
Other groups have reported6 that the initial slope of signal intensity depth profiles was determined by the extinction model, and that for larger optical depths, multiple scattering effects would cause a depth dependence governed by the reduced optical depth, rather than the optical depth.While the former result is in agreement with our finding, the nearly constant intensity level for large optical depths in intralipid (Fig. 2, left panel) differs significantly from an exponential decay governed by the reduced optical depth.
In accordance with the good localization of backscatter signal in blood even from large optical depths, estimation of Doppler frequencies from the region of blood flow showed good localization to the actual position of the vessel (Fig. 5).A significant fmding of the present work is that the average Doppler frequency for photons detected from a particular focus position was determined with an accuracy of only 3 to 4 % deviation from the mathematically expected Doppler frequency.This is remarkable, considering the fact that spectra of observed Doppler frequencies were quite broad, with standard deviations typically exceeding 300 Hz, compared to the nominal Doppler frequencies of 1210 and 1980 Hz on the flow axis for 1 5 and 25°, respectively.Apparently, the averaging of Doppler frequencies over many photons for each focus position represents such a robust estimation procedure that the resulting mean value quite accurately approximates the true frequency.
The Doppler frequency noise detected from underneath the vessel (Fig. 5) is due to the series of random Doppler shifts experienced as the photon passes through the blood flow region before and after being backscattered from static intralipid underneath the vessel.The resulting accumulated Doppler frequency for each detected photon has an expected value of zero, but a stochastic variation characterised by the standard deviation ofthe frequency distribution.
The Doppler noise level from underneath the vessel was found to be essentially independent of the incidence angle of the probe (Fig. 6), even if a larger incidence angle corresponds to higher average Doppler frequencies detected from the flow region.Furthermore, the standard deviations of Doppler noise spectra from underneath the vessel were generally greater than the standard deviations of Doppler frequencies detected from the region of blood flow (Fig. 6).This may seem surprising since the variance of the detected Doppler frequency may be expressed as a sum (over scattering locations) of contributions from individual scattering events along the paths of individual photons.Furthermore, the difference between photons backscattered from flowing blood and those backscattered from underneath the blood vessel is essentially that the latter carry no contribution to the Doppler frequency variance from the backscattering event.The reason for the smaller standard deviations of Doppler spectra from blood than from underneath the vessel must therefore be that the incoming and outgoing photon paths are more closely correlated for backscattering from blood, leading to Doppler shifts from individual scattering events that more exactly cancel each other and generate less variance when backscattering takes place within the vessel.
The conclusion from our investigations of Doppler noise is that the most precise detection of Doppler signal from a flow region relative to Doppler noise from underneath the flow region, will be achieved by a low NA detector aimed at a relatively large angle of incidence.Other effects, such as insufficient signal due to low numerical aperture, might, however, compromise this conclusion.
The present study was based on very simple modeling of the scattering properties of turbid media and does not consider inhomogeneities in refractive index or spatial correlations in the medium.Furthermore, due to computer run time limitations, simulation of full confocality of the detector was not feasible.Nevertheless, our results demonstrate the usefulness of Monte Carlo modeling for increasing the understanding of OCT and ODT, and indicate several effects from varying instrumentation parameters or medium properties on observations made by OCT and ODT.These effects are suitable for experimental verification.

Figure 3 .
Figure 3. Simulated photon tracks projected into the yz'-plane for photons detected from various depths in intralipid and blood.Detection parameters were Lc50 m, NA=O.4,O.5°: Optical depths in intralipid for panels A and C were 3 and 7 mfp units, respectively, in blood for panels B and D, 22.5 and 105 mfp units, respectively.Where visible (C, D), the lower end of the z-axis indicates the nominal focus position.

Figure 4 .
Figure 4. Mean position values for point spread functions zPSF (circles) and yPSF (triangles) for backscattering from intralipid and blood.Detection parameters were Lc=50 .tmand Lcl4 tm for intralipid and blood, respectively, with NA=O.4,8=2° for both media.The broken straight line represents the nominal focus depth position for zPSF.Error bars show I standard deviation of the PSF-distributions.For clarity, error bars are indicated as one-sided deviations for yPSF in intralipid.

3 . 2
Simulation of Doppler flow measurements Doppler flow detection was modeled according to the geometry shown in Figure1, i.e. a 100 tm diameter horizontal, cylindrical vessel was placed at an axial depth of 250 .tm in 2 % intralipid.This corresponds to an optical depth of 0.8 mfp units at the top of the vessel, and an additional depth of 15 mfp units through the blood vessel.

Figure 5 .
Figure 5. Average Doppler frequency for photons detected as a function of depth through and underneath the blood vessel modeled as shown in Figure 1 .Data are shown for probe incidence angles of I 5° and 25°.The broken curve shows the theoretically expected Doppler frequencies determined from the specified parabolic velocity profile.Detection parameters were Lc-l4 m, NA=0.4,81°.

Figure 6 .
Figure 6.Standard deviations of Doppler frequency distributions for photons detected as a function of depth through and underneath the blood vessel for NA=O.2 (triangles) and NAO.4 (circles).Other parameters as for Figure 5.
path segments that may deviate significantly from the nominal single backscattering paths.Minimally backscattered photons represent signals well-localised to the nominal focus position of the probe, resulting in narrow point spread functions with accurate mean position values.However, beyond the region of minimal scattering, detected photon intensity no longer follows the extinction model, but rather represents a constant level of diffuse photons.In this region of stray scattering the backscatter positions of detected photons no longer relate to the nominal focus position of the probe.