Open Access
11 December 2012 Effects of light scattering on optical-resolution photoacoustic microscopy
Author Affiliations +
Abstract
The imaging depth of ballistic optical imaging technologies is limited by light scattering. To study the effects of scattering on optical-resolution photoacoustic microscopy (OR-PAM), the signals were divided into target and background signals. A method to simulate the point spread function (PSF) of the PAM system considering both optical illumination and acoustic detection was proposed, then the PSF was used to calculate the contribution of each class of signal at different depths of the focal plane (zf ). How image contrast is degraded when there is a uniformly absorbing background as well as when there are small targets densely packed in the acoustic resolution cell were studied. By using the hyperboloid-focusing-based Monte Carlo method, optical focusing into a scattering medium was simulated. It was found that the lateral resolution provided by optical focusing is degraded by only 14% when zf= 1.1 transport mean free path (lt), compared with the case of no scattering. When zf= 1.7lt, the fluence at 50 μm radial distance away from the focal point is 93% of that at the focal point, which shows optical focusing is very weak at this depth. The method to simulate the PSF of PAM can be used in the future to optimize parameters so as to improve the system performance.

1.

Introduction

When using ballistic light to image in a scattering medium, such as biological tissue, light scattering attenuates the signal strength and degrades the image resolution and contrast.19 Different optical imaging modalities employ different methods to reject the background signals due to multiply scattered photons. For example, in confocal microscopy, a pinhole is used to reject multiply scattered photons and at the same time provide sectioning capability. In two-photon microscopy, the nonlinear fluorescence excitation volume is confined to the focal region. In optical coherence tomography (OCT), coherence-length gating is used to select backscattered photons with predetermined path-lengths. However, these methods are not immune to light scattering, which limits the maximum imaging depth of ballistic optical imaging technologies to about 1 mm in biological tissues.10,11 In confocal microscopy, a finite-sized pinhole accepts multiply scattered photons which degrade the sectioning capability and lateral resolution; when the size of the pinhole is too small, the signal-to-noise ratio (SNR) becomes poor.1 In two-photon microscopy, deep in-focus signal can be swamped by shallower out-of-focus fluorescence.36 In OCT, simulations have shown that image resolution is degraded with increasing depth in a highly scattering medium.8

Optical-resolution photoacoustic microscopy (OR-PAM),12 an emerging label-free optical imaging technology, is also affected by light scattering. It focuses light into biological tissue and employs an ultrasonic transducer focused on the same region to receive the photoacoustic (PA) signals generated by chromophores, showing optical absorption contrast with 100% sensitivity. Since OR-PAM employs optical focusing to provide lateral resolution, we are interested in how light scattering affects its performance.

We may describe the light scattering effects in OR-PAM mainly by the point spread function (PSF)—the image of a point object, which determines the resolution of an imaging system. Multiple scattering can severely distort the image of an embedded object, thus affecting the PSF. Therefore, both the property of the microscope and the scattering medium should be considered in quantifying the PSF.13 Usually, the PSF is a product of the illumination PSF and the detection PSF. The illumination PSF is the three-dimensional (3-D) image of a point object if the detection sensitivity is uniform in the entire space, and it is determined by the excitation light intensity distribution. The detection PSF is the 3-D image of a point object if the excitation light intensity is uniform in the entire space, and, in OR-PAM, it is determined by the sensitivity distribution of the spherically focused ultrasonic transducer.

In this paper, we present a method to simulate the PSF of OR-PAM in a scattering medium and use the PSF to study how light scattering affects the performance of OR-PAM.

2.

Method

We first defined two classes of signals in OR-PAM based on signal origin. Then, the Monte Carlo method was used to model the light propagation in biological tissue and provide the fluence distribution, which can be used to calculate the excitation of the two classes of signals. Next, the acoustic detection was simulated by taking into account the spatial-temporal response of the ultrasonic transducer. By combining the results from excitation and detection, we can obtain the PSF of the PAM system.

2.1.

Two Classes of Signals in OR-PAM

To study how light scattering degrades the strength and localization of signals, we divided the signals in OR-PAM into two classes. A Class I signal p1(t,zf) carries information about the absorption property in the target volume and hence is the desired signal (see Fig. 1):

Eq. (1)

p1(t,zf)=V1(zf)μa(z,ρ,θ)F(z,ρ,θ)p^(z,ρ,θ,t)dV,
while a Class II signal p2(t,zf) carries information from outside the target volume and obscures the Class I signal:

Eq. (2)

p2(t,zf)=V2(zf)μa(z,ρ,θ)F(z,ρ,θ)p^(z,ρ,θ,t)dV.
In Eqs. (1) and (2), t denotes the acoustic time of arrival; zf denotes the depth of the focal plane; the cylindrical coordinates z, ρ, θ represent, respectively, the depth, radial distance to the z-axis, and polar angle of a point in space; μa(z,ρ,θ) denotes the absorption coefficient; F(z,ρ,θ) denotes the optical fluence; p^(z,ρ,θ,t), whose expression will be given in Sec. 2.3, denotes the PA signal generated by a point absorber with a unit absorption coefficient and illuminated by a unit optical fluence. The region for the Class I signal, V1, corresponds to the optical resolution cell (defined by the cylinder whose base is the optical focal spot and whose height is the acoustic axial resolution), while the region for the Class II signal, V2, corresponds to the set difference of the acoustic resolution cell (defined by the cylinder whose base is the acoustic focal spot and whose height is the acoustic axial resolution) and the optical resolution cell (see Fig. 1).

Fig. 1

Definition of the Class I and the Class II signals in OR-PAM. The PA signal generated from the optical resolution cell (V1, the green region) is the Class I signal, while the PA signal generated from the set difference of the acoustic resolution cell and the optical resolution cell, i.e., the shaded region (V2), is the Class II signal.

JBO_17_12_126014_f001.png

2.2.

Simulation of Optical Fluence Distribution by the Monte Carlo Method

To simulate the illumination PSF, which reflects the light scattering effects in tissue, the Monte Carlo method is often used to calculate the optical fluence distribution. The traditional way to simulate optical focusing in the Monte Carlo method is geometric focusing, in which the initial propagation direction of a photon packet launched on the tissue surface is simply toward the geometric focus of the beam.14 This method cannot simulate the diffraction limit of the fluence distribution on the focal plane because the fluence at the focal point is always much larger than the fluence at the grid points next to the focal point.15,16 Instead, we used the hyperboloid focusing method to simulate optical focusing16,17 and we derived the formulae that can be directly used in a standard Monte Carlo simulation.11,18 In order to simulate the fluence distribution of a Gaussian beam focused into a scattering medium, we can construct the incident Gaussian beam by a set of hyperboloids of revolution of one sheet with different focal constants. A hyperboloid of one sheet is a doubly ruled surface, i.e., through each point there are two distinct lines that lie on the surface (see Fig. 2). When the position of a photon packet launched on the tissue surface is generated in the Monte Carlo simulation, the two lines passing through this point while lying on the hyperboloid can be described by their direction cosines:

Eq. (3)

Line1:ux=z^f(z^fx+y)(z^f2+1)L,uy=z^f(xz^fy)(z^f2+1)L,uz=zfL,

Eq. (4)

Line2:ux=z^f(z^fxy)(z^f2+1)L,uy=z^f(x+z^fy)(z^f2+1)L,uz=zfL,
where z^f=zf/z0 is the normalized zf, z0=(πω02)/λ is the Rayleigh range, L=(z^f2r2)/(z^f2+1)+zf2 is a normalization factor to make u=(ux,uy,uz) a unit vector, ω0=(1/π)·(λ/NA) is the beam waist of the Gaussian beam19 if the fill factor of the back aperture of the objective is 1, r is the radial position of a photon packet and x, y are the coordinates of a photon packet. Here, we set up a Cartesian coordinate system (see Fig. 2). The x-y plane is on the surface of the scattering medium, and the z axis is the normal of the surface, pointing to the scattering medium.11 The initial propagation direction of a photon packet is chosen from the two directions described by Eqs. (3) and (4) with equal probability. When there is no scattering, the trajectories of photon packets generated with the same r form a hyperboloid of revolution of one sheet. By sampling r according to the Gaussian distribution, we can generate a set of hyperboloids of one sheet. The resulting fluence distribution can perfectly reproduce the intensity distribution of a Gaussian beam when the medium is clear. In summary, we generate the position and the propagation direction of a photon packet launched on the surface of a scattering medium using the following steps:

  • 1. Generate the radial position of a photon packet by r=(ω01+z^f2)logξ1/2.

  • 2. Generate the x, y position of the photon packet from r by θ=2πξ2, x=rcos(θ), y=rsin(θ).

  • 3. Generate a random number ξ3. If 0<ξ3<0.5, the propagation direction of the photon packet is set according to Eq. (3). If 0.5ξ3<1, the propagation direction is set according to Eq. (4).

Fig. 2

Illustration of the hyperboloid focusing method. A hyperboloid of one sheet whose focal constant is 2ω0logξ1/2 (sampled from the Gaussian distribution with a 1/e2 characteristic length of ω0) is shown (in blue). The position and the two possible propagation directions of a photon packet launched on the tissue surface are shown as a red dot and two red arrows. r, x, y are the radial position and the coordinates of the photon packet. ξ1 is a random number that is uniformly distributed between 0 and 1.

JBO_17_12_126014_f002.png

In steps 1 to 3, ξi, i=1, 2, 3 is a random number that is uniformly distributed between 0 and 1. The rest of the procedure is similar to that in a standard Monte Carlo simulation. Thus, by using Eqs. (3) and (4) that we derived and following the above three steps for photon launching, the hyperboloid focusing method can be easily integrated into standard software packages such as MCML.18

In our simulations, we assumed the following tissue optical parameters:11 the scattering anisotropy g=0.9, the scattering coefficient μs=100cm1, and the absorption coefficient μa=0.1cm1. As this work targets on studying the effects of light scattering on optical focusing, we used the refractive index–matched boundary condition (i.e., the refractive indices of both the tissue and the water are 1.33) to avoid the effect of aberration caused by the mismatch between the refractive indices of the tissue and the ambient medium. The numerical aperture (NA) of the optical objective lens in air is 0.1.12 The optical wavelength was 570 nm, which is an isosbestic point for oxy- and deoxy-hemoglobin molecules. We chose grid size dr=0.1μm in the radial direction and dz=1μm in the z direction, which are more than 10 times smaller than the optical focal spot size [full width at half-maximum (FWHM) spot size=1.6μm, and the radius defined by 1/e2 of the peak value=1.36μm] and the FWHM axial resolution20 of PAM (22μm). From the optical properties of the scattering medium, we could obtain the mean free path lt=1/(μs+μa)100μm and the transport mean free path11 lt=1/[μs(1g)+μa]1mm. The illumination PSF is equal to the optical fluence distribution when there is no absorption in the scattering medium. However, we need to have some small absorption in the Monte Carlo simulation. Since the absorption length la=1/μa=10cm, which is more than 50 times larger than the maximum zf we simulated, the computed fluence distribution corresponding to each zf can be approximated as the illumination PSF for that case. The difference in fluence distribution was less than 10% when we changed μa from 0.1 to 0.001cm1.

The illumination PSF on the focal plane (illumination lateral PSF) at varied zf are shown in Fig. 3(a) and 3(b), from which we can see that the shoulder of the illumination lateral PSF rises with the increase of zf, due to scattering. When zf is close to 1lt, the FWHM of the corresponding illumination lateral PSF is not broadened much (2% when zf=0.9lt and 14% when zf=1.1lt), compared with the case of no scattering. When zf is greater than 1lt, the shoulder of the illumination lateral PSF rises very quickly with increasing zf. When zf=1.7lt, the fluence at 50 μm radial distance away from the focal point is 93% of that at the focal point, which shows that optical focusing is very weak. Due to the lack of computing power, we did not simulate the case for larger zf by directly using the Monte Carlo method. Instead, we used the diffusion theory11 to compute the illumination lateral PSF when zf=3lt, and we can see from Fig. 3(a) that no optical focusing can be discerned. The number of scattering events (Ns) for the field points on the focal plane at varied zf is shown in Fig. 3(c). Within the depth of 1lt, Ns is close to zero for field points inside the optical focal spot. When zf is greater than 1lt, Ns increases quickly. When zf=1.7lt, Ns at the focal point is greater than 40. We can also see that for a given zf, Ns for the field points outside the optical focal spot is always larger than that for the field points inside the focal spot. The on-axis fluence distributions corresponding to different zf are shown in Fig. 3(d). The fluence at the focal point decays exponentially with the increase of zf when zf is smaller than 1.3lt. The decay rate is 9.93mm1, close to μt=μs+μa=10.01mm1, which agrees with Beer’s law. Beyond 1.3lt, the fluence at the focal point decays more slowly. When zf is greater than 0.9lt, the on-axis fluence near the surface of the scattering medium becomes stronger than that on the focal plane.

Fig. 3

(a) The illumination PSF on the focal plane at varied zf. (b) A close-up of the region marked by the dash box in (a). (c) The number of scattering events for the field points on the focal plane at varied zf. (d) The on-axis fluence distributions corresponding to different zf.

JBO_17_12_126014_f003.png

2.3.

Simulation of Acoustic Detection by the Impulse Response Method

PAM normally employs a high-frequency spherically focused broadband ultrasonic transducer to detect PA signals.21 The detected pressure from a point optical absorber can be calculated as follows:22

Eq. (5)

PA(t;r)μa(r)F(r)δ(t)*SIR(t;r)*EIR(t),
where * denotes temporal convolution, r is the position of the PA point source, EIR(t) is the electro-mechanical impulse response of the transducer, and SIR(t;r)=S0{δ[t(|rr0|/c)]}/(2π|rr0|)dS0 (where S0 denotes the aperture of the transducer, and r0 denotes the position of the point on the aperture) is the spatial impulse response of the transducer, which describes the diffraction of sound by the aperture of the transducer.23 It is important to note that SIR here is computed for an ultrasonic transducer that is focused by curving the piezoelectric element directly instead of by attaching an acoustic lens. In the latter case, we can compute the SIR of an equivalent curved transducer whose focal length is the same as that of the acoustic lens. δ'(t) is the time derivative of the delta function and is the PA signal generated by a unit (i.e., μaF=1) point absorber.11 δ(t)*SIR(t;r)*EIR(t) is the expression for p^(z,ρ,t) in Eqs. (1) and (2). δ(t)*EIR(t) can be determined experimentally by measuring the PA signal of a point absorber at the focal point of the ultrasonic transducer, where SIR is a delta function of time.24 For simulation purposes, we generate the EIR numerically by Eq. (6):

Eq. (6)

EIR(t)=sin(2πf0t)·exp(t22a2),
where t is time, a=(2ln2/π)·(1/FWHMf), f0 is the center frequency of the transducer, and FWHMf is the 6dB bandwidth (one way for receiving only). We chose f0=50MHz, FWHMf=60MHz in our simulations. The generated EIR is shown in Fig. 4(a) and its spectrum is shown in Fig. 4(b), from which we can see that the shape of the spectrum is close to Gaussian. We can calculate the waveform of the pulse-echo response of the transducer by convolving EIR with itself. The calculated waveform and its spectrum (data not shown) resemble well the result in the test report of the transducer (Olympus NTD V214). The duration of EIR as measured between the outer zero-crossing points of the two major lobes is about 20 ns. In Fig. 4(a), δ(t)*EIR(t) is also shown.

Fig. 4

(a) EIR(t) generated by Eq. (6) and (b) its spectrum. In (a),δ(t)*EIR(t) is also shown.

JBO_17_12_126014_f004.png

We calculated SIR(t;r) of the field points near or on the focal plane in the frequency domain as follows:25

Eq. (7a)

SIR(ω;r)=eik|rr0|unA02π|rr0|G(z),

Eq. (7b)

G(z)=(2/z)n=0(1)n(d0/a0)2nJ2n+1(z),

Eq. (7c)

z=[1i/(k|rr0|)]ka0sinθka0sinθ,
where k=ω/c is the wavenumber, c is the speed of sound in the medium, a0 is the radius of the transducer (see Fig. 5), d0 is the depth of the concave surface, A0 is the area of the aperture, un is the normal velocity of the particle on the aperture, and θ is the angle from the axis of the transducer to the line connecting the apex of the aperture to a field point. J2n+1(z) is the Bessel function. The duration of SIR of a field point near the focal plane is Δt=(r2r1)/c, where r1 and r2 denote the distances from a field point to the closest and furthest edge (E1 and E2, respectively) of the transducer24 (see Fig. 5). The NA of our transducer is 0.58 and its focal length is 5.17 mm, from which we can get a0=3mm and d0=0.96mm. Though the author assumed a low NA in the derivation25 of Eq. (7), Coulouvrat26 has shown that the formula also applies to the case of a high NA, which enables our use of Eq. (7) in this work. The 500 MHz sampling rate of our OR-PAM system, providing proper sampling of PA signals up to 250 MHz according to the Nyquist sampling theorem, determines the maximum frequency that we need to simulate SIR. In addition, the maximum duration of SIR of all simulated field points determines the sampling interval in the frequency domain (2.2 MHz was used in our simulations). To make the time domain signal smoother in Figs. 6(b) and 7(a), we zero-padded SIR(ω;r) to 500 MHz before taking the inverse discrete Fourier transform to interpolate the signal in the time domain.

Fig. 5

Schematic of a spherically focused transducer. a0 is the radius of the transducer, d0 is the depth of the concave surface, f is the focal length of the transducer, θ is the angle from the axis of the transducer to the line connecting the apex of the aperture to a field point P, r1 and r2 are the distances from P to the closest and furthest edge (E1 and E2, respectively) of the transducer.

JBO_17_12_126014_f005.png

Fig. 6

The calculated SIR of four PA point sources on the focal plane (a) in the frequency domain and (b) in the time-domain. ρ is the radial distance between the point source and the focal point. The solid arrow in (b) denotes the delta function. In (a), the spectrum of EIR is also plotted for comparison.

JBO_17_12_126014_f006.png

Fig. 7

(a) Simulated detected PA signals of four point sources on the focal plane. ρ is the radial distance between the point source and the focal point. (b) The detection PSF on the focal plane.

JBO_17_12_126014_f007.png

The magnitude of the calculated SIR(ω;r) for four PA point sources on the focal plane is shown in Fig. 6(a), and the SIR in the time domain is shown in Fig. 6(b). It can be seen that the further the PA point source is from the focal point, the narrower the bandwidth of its SIR, and the longer the duration of SIR in the time domain. When a PA source on the focal plane is 40 μm away from the focal point, the spectrum of its SIR has little overlap with the spectrum of EIR. Moreover, in the time domain, the corresponding SIR is broadened and its duration is comparable to the duration of EIR. Because of the convolution in Eq. (5), this will cause an obvious distortion of the waveform of the detected PA signal, which can be seen from Fig. 7(a).

By using the frequency domain formula of Eq. (5) and setting μa(r)F(r)=1, the detected PA signals generated by the same four point sources as in Fig. 6(a) were calculated and the converted time-domain signals are shown in Fig. 7(a). By plotting the amplitude of the PA signal of each point source at the time when the signal generated from the focal point reaches its maximum [denoted by tm in Fig. 7(a)], the detection PSF on the focal plane (detection lateral PSF) is obtained and shown in Fig. 7(b). Unlike the sensitivity curve of a transducer working in continuous-wave (CW) mode, there are no sidelobes in the detection lateral PSF for broadband detection. The radius of the acoustic focal spot (defined as 1/e2 of the peak value) is 24.5 μm, and the FWHM spot size is 29.2 μm.

3.

Results

3.1.

Point Spread Function, Line Spread Function, and Edge Spread Function of PAM

Since the depth of focus of the acoustic focal zone (220μm) is more than ten times larger than the axial resolution of PAM (22μm), the detection PSF at each depth in the acoustic resolution cell is nearly the same as that on the focal plane (data not shown). By multiplying the illumination PSF and the detection PSF, the PSF of PAM is obtained. The PSFs on the focal plane (lateral PSF) at varied zf are shown in Fig. 8(a). The line spread function (LSF) and the edge spread function (ESF) were calculated from the PSF11 and are shown in Fig. 8(b) and 8(c). The LSF and the ESF represent the 1-D image of a line and an edge, respectively.

Fig. 8

(a) The lateral PSF, (b) LSF, and (c) ESF of PAM at varied zf. The horizontal axis is perpendicular to the line and the edge.

JBO_17_12_126014_f008.png

From Fig. 8(a), we can see that for zf smaller than 1.3lt, optical focusing is finer than acoustic focusing. The PAM working within this depth range is called the optical-resolution PAM.12,27 For zf greater than 1.3lt, optical focusing becomes very weak and the resolution is mainly determined by acoustic focusing. The PAM working within this depth range is called the acoustic-resolution PAM (AR-PAM).27,28 From Fig. 8, we can see that the lateral PSF, LSF, and ESF are degraded by a base with increasing zf, due to scattering. Moreover, for a given zf, the LSF is always broader than the PSF. This is because no Class II signal exists when only a point target is present; however, for a line target, the Class II signal becomes nonzero according to Eq. (2). We can also see in Fig. 8(c) that, for zf greater than 1.1lt, the ESFs are visually identical. However, in Fig. 8(b), we can still distinguish the LSF for zf=1.1lt from that for zf=3lt. In this sense, we may say that the contrast of the ESF is worse than that of the LSF for a given zf. This is due to the fact that the Class II signal becomes stronger when the object changes from a line to an edge, which can be seen from Eq. (2).

3.2.

Contrast Degradation When Imaging a Target in a Uniform Background

We investigated the effects of light scattering when imaging a target in a uniformly absorbing background. Let us consider the task of using PAM to image cell nuclei29,30 in the background of protein [see Fig. 9(a)]. When light is focused on the target (DNA/RNA), the measured signal can be written as

Eq. (8)

MT(zf)=k0[μa1PA1(zf)+μa2PA2(zf)],
where μa1 and μa2 are the absorption coefficients of the target and the background, k0 is a constant factor, and PA1 and PA2 are the peak amplitudes of the Class I and the Class II signals. Hereafter, the Class I and the Class II signals refer to their peak amplitudes.

Fig. 9

(a) Schematic of imaging a target in a uniformly absorbing background. (b) The Class I and the Class II signals at varied zf and their ratios.

JBO_17_12_126014_f009.png

When light is focused on the background (protein), the measured signal can be approximated as follows by ignoring the contribution from the target:

Eq. (9)

MB(zf)=k0[μa2PA1(zf)+μa2PA2(zf)].
We define the ratio of the Class I signal to the Class II signal Rp(zf)=PA1(zf)/PA2(zf) and the ratio of the absorption coefficient of the target to that of the background Ra=μa1/μa2. Then, the measurement ratio RM(zf)=MT(zf)/MB(zf) can be expressed as

Eq. (10)

RM(zf)=RaRa1Rp(zf)+1
=Rp(zf)·RaRp(zf)+1+1Rp(zf)+1

Eq. (11)

{Ra,whenRp(zf)1Rp(zf)·Ra+1,whenRp(zf)1.
From Eq. (10), we can see that the measurement ratio is always smaller than the ratio of absorption coefficients, which can be regarded as the true contrast. Equation (11) shows that when the Class II signal is much weaker than the Class I signal, the measurement ratio can almost recover the ratio of absorption coefficients, thus providing nearly true contrast. However, when the Class II signal is much stronger than the Class I signal, the measurement ratio cannot recover the ratio of absorption coefficients and will degrade the contrast. In fact, the contrast to noise ratio11 (CNR) is related to RM as follows:

Eq. (12)

CNR(zf)=MT(zf)MB(zf)MB(zf)σN=RM(zf)1σN

Eq. (13)

{RaσN,whenRp(zf)1andRa1RaσN·Rp(zf),whenRp(zf)1.
In Eqs. (12) and (13), σN is the standard deviation of the background intensity. From Eqs. (10) and (12), we can see that the CNR decreases with the increase of the proportion of the Class II signal in the total signal. If Ra1, which is the case for thymus DNA and glutamate dehydrogenase (Ra35 at 250 nm29), Eq. (13) shows that when the Class II signal is much stronger than the Class I signal, the CNR will decrease by a factor of 1/Rp(zf)=PA2(zf)/PA1(zf), compared with the CNR in the absence of scattering (Rp=, CNR=Ra/σN).

The Class I and the Class II signals at varied zf and their ratios are shown in Fig. 9(b). It can be seen that the Class II signal is stronger than the Class I signal for zf greater than 0.9lt (0.09lt). When zf is smaller than 1.1lt, the signal decay rate of the Class I signal is 9.87mm1, close to μt=10.01mm1, which agrees with Beer’s law. In contrast, the signal decay rate of the Class II signal is 5.28mm1, which means the Class II signal decays much more slowly than the Class I signal. Since the CNR is closely related to Rp, we can get an idea of how it is degraded with the increase of zf by Eq. (13). For example, Rp is 1/19 when zf=0.9lt and 1/130 when zf=1.7lt which means a decrease of 19 and 130 times in CNR compared with that when there is no scattering. It should be noted that in this model we assumed the target size is comparable to the lateral resolution in the absence of scattering. When the target size is different, the model will be more complicated and Eqs. (1) and (2) should be modified.

3.3.

Contrast Degradation When Imaging Small Targets Densely Packed in the Acoustic Resolution Cell

To study the effects of light scattering on imaging small targets that are densely packed in the acoustic resolution cell, such as capillaries in a capillary bed, an array of evenly spaced capillaries were used as a simple model for the capillary bed [see Fig. 10(a)]. To change the packing density of the target, we varied the distance between capillaries from 8 to 113 μm, while keeping the diameter of the capillary unchanged (8 μm). Using the lateral PSF, we can simulate the image of a 2-D object O(x,y) in a scattering medium by convolution:13

Eq. (14)

I(x,y)=+O(x,y)PSF(xx,yy)dxdy.
Figure 10(b) and 10(c) show the images of capillaries when the distance between two capillaries is 8 and 113 μm, from which we can calculate the modulation or contrast31 C=(ImaxImin)/(Imax+Imin). The contrasts corresponding to different distances between two capillaries (8, 24, 40, 56, and 113 μm) and different zf are shown in Fig. 10(d). It can be seen that the contrast is degraded with increasing zf and target density. Different from the scenario in Fig. 9(a) where the Class II signal comes from everywhere in the region for the Class II signal (V2 in Fig. 1), here, the Class II signal is generated from PA sources that are densely packed in V2. When the packing density of the acoustic resolution cell is increased, more Class II signals will be generated according to Eq. (2), which explains why the contrast decreases with the increase of packing density. It should be noted that the scenario considered here is a simplified model, because in reality, the capillaries in a capillary bed may not lie on the same plane perpendicular to the optical axis and the distance between two capillaries may not be the same for different pairs of capillaries.

Fig. 10

(a) Schematic of imaging an array of evenly spaced capillaries. (b) The image of capillaries when the distance between two capillaries is the width of the capillary (8 μm) and (c) 113 μm. (d) The contrasts corresponding to different zf and different distances between two capillaries.

JBO_17_12_126014_f010.png

4.

Discussion

Optical focusing in a scattering medium is highly desired in optical imaging and sensing, photodynamic therapy and manipulation.32 The transport mean free path, lt, indicates the mean propagation distance that it takes for photons to lose memory of the initial propagation direction they had before entering the scattering medium.10 Thus, optical focusing beyond 1lt is generally regarded as infeasible. However, the lateral resolution when zf=1lt and how it decays with depth have not been studied sufficiently, either theoretically or experimentally. In this paper, we found that when zf is close to 1lt, the FWHM of the corresponding illumination lateral PSF is not broadened much (2% when zf=0.9lt and 14% when zf=1.1lt), compared with the case of no scattering. This seems somewhat surprising, yet agrees with the simulation results by Hayakawa et al.33 When zf is greater than 1lt, the shoulder of the illumination lateral PSF rises very quickly with increasing zf. When zf=1.7lt, the fluence at 50 μm radial distance away from the focal point is 93% of that at the focal point, which shows that optical focusing is very weak at this depth. Experiments to validate these results are expected to be done in future work. It is worth mentioning that it has recently been proposed to phase conjugate the ultrasound-modulated photons to form an optical focus inside a thick turbid medium.32 With this so-called time-reversed ultrasonically encoded (TRUE) optical focusing technique, researchers have thus far successfully demonstrated dynamic optical focusing in tissue-mimicking phantoms and ex vivo tissue samples with thicknesses up to 12 to 16lt34,35 and 10lt,36 respectively, far exceeding the 1lt limit of passive optical focusing. Through digital-optical-phase-conjugation-based TRUE focusing, fluorescence imaging beyond the ballistic regime has subsequently been realized.37,38

From the on-axis fluence distribution in Fig. 3(d), we can see that when zf is greater than 0.9lt, the fluence near the surface of the scattering medium becomes stronger than that on the focal plane, which degrades the image contrast of two-photon microscopy and is known as the fundamental limit of the maximum imaging depth of this image modality and other nonlinear optical microscopy.4,5 In contrast, because of the time-resolved acoustic detection, the contrast of a PAM image is not affected by the strong near-surface fluence but is mainly degraded by the Class II signal. In comparison, in two-photon microscopy, because its illumination lateral PSF is that of PAM squared, the contrast is less affected by the Class II signal around the focal plane. Considering the above analysis, it might be a good idea to develop two-photon PAM39 to improve the image contrast of PAM. The use of longer wavelength may also increase the maximum imaging depth. However, the low efficiency of two-photon absorption and how to separate the two-photon PA signal from the single-photon PA signal are challenging.

In this work, we simulated focusing light into a semi-infinite scattering medium, since we intended to study the maximum imaging depth of OR-PAM. In practice, OR-PAM may be used to image thin tissues, such as the mouse ear,40 cells,41 and zebrafish larva.42 In these cases, the bottom boundary of the thin tissue may have an effect on the light distribution. However, the effect on optical focusing is minor because most reflected/scattered light from the bottom boundary cannot reach the optical focus.

It is important to note that the Monte Carlo method only describes the transport of energy, therefore it is incapable of modeling coherent phenomena. This limitation also applies to the hyperboloid-focusing-based Monte Carlo method. Although an imperfect method to simulate optical focusing in a scattering medium, it has been validated by some models and experiments. In the original paper that described this method,16 the heterodyne efficiency factors and the transverse intensity distribution simulated by the hyperboloid-focusing-based Monte Carlo method agree with those simulated by the extended Huygens-Fresnel (EHF) model,43 which is a widely acknowledged model in OCT and has been validated experimentally. Good agreement between the heterodyne efficiency factors obtained by the Monte Carlo simulation and by experiments was also reported.44 Moreover, two-photon fluorescence signals were simulated by the hyperboloid-focusing-based Monte Carlo method, and the result agreed well with the experiments.3

The hyperboloid focusing method used in the Monte Carlo simulation provides a way to simulate optical focusing. However, when the grid size is small (e.g., dr=0.1μm and dz=1μm were used in this work), in order to obtain an acceptable statistical error in the result, 2×1010 photon packets with a cut-off weight of 104 were used to simulate the fluence distribution when zf=0.9lt, and 4×1012 photon packets were used when zf=1.7lt. Thus, graphics processing units (GPUs) are highly recommended to accelerate the Monte Carlo simulation.45

5.

Conclusions

In this work, we employed the hyperboloid focusing method in the Monte Carlo simulation to simulate optical focusing in a scattering medium. We also used the impulse response method to simulate the detection of PA signal by a spherically focused broadband high-NA ultrasonic transducer. By combining the results from excitation and detection, we simulated the PSF of the PAM system. The PSF can be used in the future to optimize parameters so as to improve the system performance. For example, we can study the effects of the transducer bandwidth on the lateral resolution of AR-PAM. To study the effects of light scattering, we divided the signals in OR-PAM into two classes and used the PSF to calculate the contribution of each class of signal to the total signal at varied zf. We studied quantitatively how image contrast is degraded when there is a uniformly absorbing background as well as when there are small targets densely packed in the acoustic resolution cell. Moreover, we simulated optical focusing into a scattering medium, especially when the range of zf was from 1lt to 1.7lt. It was found that the lateral resolution provided by optical focusing is degraded by only 14% when zf is increased to 1.1lt, compared with the case of no scattering. However, the lateral resolution was degraded very quickly with further increase of zf. When zf=1.7lt, the fluence at 50 μm radial distance away from the focal point is 93% of that at the focal point, which shows that optical focusing is very weak at this depth. Experimental validation of these simulation results is expected in future work.

Acknowledgments

We would like to thank Zijian Guo, Kun Wang, and Lidai Wang for helpful discussions. We would also like to thank the Cloud cluster of Washington University for computing support. This work was sponsored in part by National Institutes of Health grants R01 EB008085, R01 EB000712, R01 CA134539, R01 CA157277, R01 CA159959, and U54 CA 136398. L.W. has a financial interest in Microphotoacoustics, Inc. and Endra, Inc., which, however, did not support this work.

References

1. 

J. M. SchmittA. KnuttelM. Yadlowsky, “Confocal microscopy in turbid media,” J. Opt. Soc. Am. A., 11 (8), 2226 –2235 (1994). http://dx.doi.org/10.1364/JOSAA.11.002226 JOAOD6 1084-7529 Google Scholar

2. 

J. P. YingF. LiuR. R. Alfano, “Spatial distribution of two-photon-excited fluorescence in scattering media,” Appl. Opt., 38 (10), 2151 (1999). http://dx.doi.org/10.1364/AO.38.002151 APOPAI 0003-6935 Google Scholar

3. 

N. J. Durret al., “Maximum imaging depth of two-photon autofluorescence microscopy in epithelial tissues,” J. Biomed. Opt., 16 (2), 026008 (2011). http://dx.doi.org/10.1117/1.3548646 JBOPFO 1083-3668 Google Scholar

4. 

P. TheerW. Denk, “On the fundamental imaging-depth limit in two-photon microscopy,” J. Opt. Soc. Am. A, 23 (12), 3139 –3149 (2006). http://dx.doi.org/10.1364/JOSAA.23.003139 JOAOD6 0740-3232 Google Scholar

5. 

P. Theer, “On the fundamental imaging-depth limit in two-photon microscopy,” University of Heidelberg, (2004). Google Scholar

6. 

P. TheerM. T. HasanW. Denk, “Two-photon imaging to a depth of 1000 μm in living brains by use of a TiAl2O3 regenerative amplifier,” Opt. Lett., 28 (12), 1022 –1024 (2003). http://dx.doi.org/10.1364/OL.28.001022 OPLEDP 0146-9592 Google Scholar

7. 

D. J. Smithieset al., “Signal attenuation and localization in optical coherence tomography studied by Monte Carlo simulation,” Phys. Med. Biol., 43 (10), 3025 –3044 (1998). http://dx.doi.org/10.1088/0031-9155/43/10/024 PHMBA7 0031-9155 Google Scholar

8. 

R. Wang, “Signal degradation by coherence tomography multiple scattering in optical of dense tissue: a Monte Carlo study towards optical clearing of biotissues,” Phys. Med. Biol., 47 (13), 2281 –2299 (2002). http://dx.doi.org/10.1088/0031-9155/47/13/307 PHMBA7 0031-9155 Google Scholar

9. 

G. YaoL. V. Wang, “Monte Carlo simulation of an optical coherence tomography signal in homogeneous turbid media,” Phys. Med. Biol., 44 (9), 2307 –2320 (1999). http://dx.doi.org/10.1088/0031-9155/44/9/316 PHMBA7 0031-9155 Google Scholar

10. 

V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods, 7 (8), 603 –614 (2010). http://dx.doi.org/10.1038/nmeth.1483 1548-7091 Google Scholar

11. 

L. V. WangH. Wu, Biomedical Optics : Principles and Imaging, Wiley-Interscience, Hoboken, NJ (2007). Google Scholar

12. 

K. Maslovet al., “Optical-resolution photoacoustic microscopy for in vivo imaging of single capillaries,” Opt. Lett., 33 (9), 929 –931 (2008). http://dx.doi.org/10.1364/OL.33.000929 OPLEDP 0146-9592 Google Scholar

13. 

X. S. GanM. Gu, “Effective point-spread function for fast image modeling and processing in microscopic imaging through turbid media,” Opt. Lett., 24 (11), 741 –743 (1999). http://dx.doi.org/10.1364/OL.24.000741 OPLEDP 0146-9592 Google Scholar

14. 

L. V. WangG. Liang, “Absorption distribution of an optical beam focused into a turbid medium,” Appl. Opt., 38 (22), 4951 –4958 (1999). http://dx.doi.org/10.1364/AO.38.004951 APOPAI 0003-6935 Google Scholar

15. 

Y. Liuet al., “Effect of light scattering on optical-resolution photoacoustic microscopy,” Proc. SPIE, 8223 82233S (2012). http://dx.doi.org/10.1117/12.909373 APOPAI 0003-6935 Google Scholar

16. 

A. Tychoet al., “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt., 41 (31), 6676 –6691 (2002). http://dx.doi.org/10.1364/AO.41.006676 APOPAI 0003-6935 Google Scholar

17. 

F. Zhanget al., “Monte Carlo method for simulating optical coherence tomography signal in homogeneous turbid media,” Proc. SPIE, 7022 702213 (2007). http://dx.doi.org/10.1117/12.804101 PSISDG 0277-786X Google Scholar

18. 

L. V. WangS. L. JacquesL. Q. Zheng, “MCML—Monte-Carlo modeling of light transport in multilayered tissues,” Comput. Meth. Prog. Bio., 47 (2), 131 –146 (1995). http://dx.doi.org/10.1016/0169-2607(95)01640-F CMPBEK 0169-2607 Google Scholar

19. 

B. E. A. SalehM. C. Teich, Fundamentals of Photonics, Wiley-Interscience, Hoboken, NJ (2007). Google Scholar

20. 

C. Zhanget al., “In vivo photoacoustic microscopy with 7.6-μm axial resolution using a commercial 125-MHz ultrasonic transducer,” J. Biomed. Opt., 17 (11), 116016 (2012). http://dx.doi.org/10.1117/1.JBO.17.11.116016 JBOPFO 1083-3668 Google Scholar

21. 

H. F. Zhanget al., “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nat. Biotechnol., 24 (7), 848 –851 (2006). http://dx.doi.org/10.1038/nbt1220 NABIF9 1087-0156 Google Scholar

22. 

K. Wanget al., “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol., 57 (17), 5399 –5423 (2012). http://dx.doi.org/10.1088/0031-9155/57/17/5399 PHMBA7 0031-9155 Google Scholar

23. 

P. R. Stepanishen, “Transient radiation from pistons in an infinite planar baffle,” J. Acoust. Soc. Am., 49 (5), 1629 (1971). http://dx.doi.org/10.1121/1.1912541 JASMAN 0001-4966 Google Scholar

24. 

M. ArditiF. S. FosterJ. W. Hunt, “Transient fields of concave annular arrays,” Ultrason. Imag., 3 (1), 37 –61 (1981). http://dx.doi.org/10.1016/0161-7346(81)90081-X ULIMD4 0161-7346 Google Scholar

25. 

H. T. O’Neil, “Theory of focusing radiators,” J. Acoust. Soc. Am., 21 (5), 516 –526 (1949). http://dx.doi.org/10.1121/1.1906542 JASMAN 0001-4966 Google Scholar

26. 

F. Coulouvrat, “Continuous field radiated by a geometrically focused transducer—numerical investigation and comparison with an approximate model,” J. Acoust. Soc. Am., 94 (3), 1663 –1675 (1993). http://dx.doi.org/10.1121/1.408139 JASMAN 0001-4966 Google Scholar

27. 

L. V. WangS. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science, 335 (6075), 1458 –1462 (2012). http://dx.doi.org/10.1126/science.1216210 SCIEAS 0036-8075 Google Scholar

28. 

K. MaslovG. StoicaL. V. Wang, “In vivo dark-field reflection-mode photoacoustic microscopy,” Opt. Lett., 30 (6), 625 –627 (2005). http://dx.doi.org/10.1364/OL.30.000625 OPLEDP 0146-9592 Google Scholar

29. 

D. K. Yaoet al., “Optimal ultraviolet wavelength for in vivo photoacoustic imaging of cell nuclei,” J. Biomed. Opt., 17 (5), 056004 (2012). http://dx.doi.org/10.1117/1.JBO.17.5.056004 JBOPFO 1083-3668 Google Scholar

30. 

D. K. Yaoet al., “In vivo label-free photoacoustic microscopy of cell nuclei by excitation of DNA and RNA,” Opt. Lett., 35 (24), 4139 –4141 (2010). http://dx.doi.org/10.1364/OL.35.004139 OPLEDP 0146-9592 Google Scholar

31. 

E. Hecht, Optics, Addison-Wesley, Reading, MA (2002). Google Scholar

32. 

X. XuH. LiuL. V. Wang, “Time-reversed ultrasonically encoded optical focusing into scattering media,” Nat. Photon., 5 (3), 154 –157 (2011). http://dx.doi.org/10.1038/nphoton.2010.306 1749-4885 Google Scholar

33. 

C. K. Hayakawaet al., “Amplitude and phase of tightly focused laser beams in turbid media,” Phys. Rev. Lett., 103 (4), 043903 (2009). http://dx.doi.org/10.1103/PhysRevLett.103.043903 PRLTAO 0031-9007 Google Scholar

34. 

P. Laiet al., “Reflection-mode time-reversed ultrasonically encoded optical focusing into turbid media,” J. Biomed. Opt., 16 (8), 080505 (2011). http://dx.doi.org/10.1117/1.3609001 JBOPFO 1083-3668 Google Scholar

35. 

Y. Suzukiet al., “Energy enhancement in time-reversed ultrasonically encoded optical focusing using a photorefractive polymer,” J. Biomed. Opt., 17 (8), 080507 (2012). http://dx.doi.org/10.1117/1.JBO.17.8.080507 JBOPFO 1083-3668 Google Scholar

36. 

P. Laiet al., “Time-reversed ultrasonically encoded optical focusing in biological tissue,” J. Biomed. Opt., 17 (3), 030506 (2012). http://dx.doi.org/10.1117/1.JBO.17.3.030506 JBOPFO 1083-3668 Google Scholar

37. 

K. SiR. FiolkaM. Cui, “Fluorescence imaging beyond the ballistic regime by ultrasound-pulse-guided digital phase conjugation,” Nat. Photon., 6 (10), 657 –661 (2012). http://dx.doi.org/10.1038/nphoton.2012.205 1749-4885 Google Scholar

38. 

Y. Wanget al., “Deep-tissue focal fluorescence imaging with digitally time-reversed ultrasound-encoded light,” Nat. Commun., 3 928 (2012). http://dx.doi.org/10.1038/ncomms1925 Google Scholar

39. 

Y. YamaokaM. NambuT. Takamatsu, “Fine depth resolution of two-photon absorption-induced photoacoustic microscopy using low-frequency bandpass filtering,” Opt Express, 19 (14), 13365 –13377 (2011). http://dx.doi.org/10.1364/OE.19.013365 OPEXFF 1094-4087 Google Scholar

40. 

J. J. Yaoet al., “Label-free oxygen-metabolic photoacoustic microscopy in vivo,” J. Biomed. Opt., 16 (7), 076003 (2011). http://dx.doi.org/10.1117/1.3594786 JBOPFO 1083-3668 Google Scholar

41. 

C. ZhangK. MaslovL. V. Wang, “Subwavelength-resolution label-free photoacoustic microscopy of optical absorption in vivo,” Opt. Lett., 35 (19), 3195 –3197 (2010). http://dx.doi.org/10.1364/OL.35.003195 OPLEDP 0146-9592 Google Scholar

42. 

S. Yeet al., “Label-free imaging of zebrafish larvae in vivo by photoacoustic microscopy,” Biomed. Opt. Express, 3 (2), 360 –365 (2012). http://dx.doi.org/10.1364/BOE.3.000360 BOEICL 2156-7085 Google Scholar

43. 

L. ThraneH. T. YuraP. E. Andersen, “Analysis of optical coherence tomography systems based on the extended Huygens-Fresnel principle,” J. Opt. Soc. Am. A, 17 (3), 484 –490 (2000). http://dx.doi.org/10.1364/JOSAA.17.000484 JOAOD6 0740-3232 Google Scholar

44. 

L. Thraneet al., “Extraction of optical scattering parameters and attenuation compensation in optical coherence tomography images of multilayered tissue structures,” Opt. Lett., 29 (14), 1641 –1643 (2004). http://dx.doi.org/10.1364/OL.29.001641 OPLEDP 0146-9592 Google Scholar

45. 

E. Alerstamet al., “Next-generation acceleration and code optimization for light transport in turbid media using GPUs,” Biomed. Opt. Express, 1 (2), 658 –675 (2010). http://dx.doi.org/10.1364/BOE.1.000658 BOEICL 2156-7085 Google Scholar
© 2012 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2012/$25.00 © 2012 SPIE
Yan Liu, Chi Zhang, and Lihong V. Wang "Effects of light scattering on optical-resolution photoacoustic microscopy," Journal of Biomedical Optics 17(12), 126014 (11 December 2012). https://doi.org/10.1117/1.JBO.17.12.126014
Published: 11 December 2012
Lens.org Logo
CITATIONS
Cited by 76 scholarly publications.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Point spread functions

Monte Carlo methods

Scattering

Light scattering

Tissue optics

Acoustics

Transducers

Back to Top