2.1.

2PAM Contrast

In the absence of fluorescence saturation, photobleaching, and sample movement, the effective fluorophore concentration is independent of time. Assuming the concentration of fluorophores within V _s is a constant, C _s, and that the out-of-focus fluorophore concentration is diffuse enough to be approximated by the average fluorophore concentration, C _b, we can remove the concentration terms from the integral in Eq. 1. We define the staining inhomogeneity, χ, as the ratio of the fluorophore concentrations, and R as the ratio of integrated intensities:

We can then express the ratio of signal to background fluorescence as the product of χ and R:

Though our definition for F _s/F _b is the same as previous reports, our definition of χ and R differ from that described previously.¹⁸ Our definition of the staining inhomogeneity relies on the absolute concentrations of fluorophores inside and outside a signal volume and does not depend on the perifocal intensity distribution (i.e., the numerical aperture). Additionally, our R value is defined by the volume of a feature of interest rather than the volume of the focal spot. Our approach has the distinction that the χ value is an intrinsic, system-independent property of the sample that scales linearly with changes in background or signal fluorophore concentration.

In this paper, we explore the gradual decay of image contrast, Q, versus depth due to out-of-focus fluorescence. We define the contrast in terms of the measured signal (MS) when the focal spot concentrically overlaps with V _s and the measured background (MB) when the focal spot is in a volume at the same imaging depth that produces a relatively weak signal (Fig. 1), such that:

Fig. 1

Diagram of the parameters used in quantifying the contrast in 2PAM. The signal and background volumes (V _s and V _b, respectively) and signal and background fluorophore concentrations (C _s and C _b, respectively) are indicated. MS is defined as the fluorescence detected when the geometric focal point of the excitation light is at the center of V _s. MB is defined as the fluorescence detected when the excitation is focused at the same imaging depth as MS, but outside of the feature of interest.

2.2.

2PAM System

We built a two-photon microscope specifically designed for deep tissue autofluorescence imaging. An 80-MHz repetition rate, Ti:Saphire laser oscillator (Spectra Physics, Mai Tai) generated the excitation pulses with energies of up to 10 nJ. We used 760-nm excitation wavelength because it yielded the brightest autofluorescence signal with minimal excitation light bleed-through in our system. Two half-waveplate/polarizing beam cube pairs attenuated the excitation beam and a pair of galvanometer scanning mirrors (Cambridge Technologies, 6215h) scanned the laser into an upright microscope. We used a high numerical aperture (NA), water dipping objective with a large field of view and a 2-mm working distance (0.95/20x Olympus XLUMPFL). The large field of view is especially important for high collection efficiency at large imaging depths.²⁷ The physical radii of the back aperture, r _BA, and front aperture, r _FA, of this objective are 8.5 and 2 mm, respectively. The measured beam 1/e² radius of the excitation beam, w _BA, was 6 mm at the objective back aperture, giving a fill factor of 0.7. The transmission of 760-nm light was measured to be 70% by measuring the laser power before and after the objective with a power meter.

We optimized the collection optics for maximum collection of emission photons from the diffusive regime using Zemax-EE (2009 Version), similar to a previous approach.²⁸ Excess stray and scattered excitation light were filtered through two low pass filters (3-mm thick BG-39, Schott Glass and FF01–750, Semrock). We used a high-sensitivity, cooled, GaAsP photomultiplier tube (PMT) module (Hamamatsu, H7422–40) as our detector. The current from the PMT was sent through a 1-MHz bandwidth transimpedance amplifier (Stanford Research Systems SR570) and sampled at 1 MHz with 12-bit digitization on a data acquisition card (National Instruments PCI-6115). We used a modified version of MPScan software²⁹ to control the laser power and scanning mirrors, and acquire 512×512 pixel images at 3 frames/s. All images presented and analyzed in this paper are the average of 8 frames to minimize other sources of noise. In biopsy and phantom measurements, we increased the excitation power manually with depth to maintain an approximately constant level of collected signal at levels just below the saturation limit of the system.

We characterized the spatial intensity distribution of the excitation light at the focal volume by measuring the intensity-squared point spread function (IPSF²) of our system with 100-nm diameter fluorescent beads embedded in agar. A representative IPSF² measurement is shown in Fig. 2a. We found lateral and axial FWHMs of r _FWHM = 460±60 nm and z _FWHM = 1760±130 nm (mean ± standard deviation), respectively, by fitting Gaussian functions to the profiles drawn through the centroids of 20 beads. These values are approximately 50% larger than theoretical values expected with a diffraction-limited spot from a 0.95 NA water dipping objective,³⁰ and closer to what would be expected from diffraction-limited focusing from an NA of 0.75. The difference between measured and diffraction-limited IPSF² is too large to be accounted for by our underfilling of the back aperture, ^{31, 32} and is likely due to lens aberrations at NIR wavelengths. Similarly large IPSF² have been previously reported from this objective. ^{27, 33, 34} We found that the IPSF² was independent of imaging depth through the full working distance of the lens (2 mm) in a sample of 100-nm fluorescent beads embedded in an agar gel. The constant spot size indicates that specimen-induced aberrations negligibly affect the shape of the intensity distribution in the perifocal volume in agar phantoms, which is consistent with other studies. ^{15, 16, 35}

Fig. 2

A representative lateral, r, and axial, z, ISPF² from a 100-nm fluorescent bead (a), indicating average lateral and axial FWHMs of 460 and 1760 nm, respectively. The inset shows the two-photon image of a typical bead in the rz-plane. Autocorrelation measurements with a 1 in. singlet lens (b) and with a 20x/0.95 objective lens (c) indicated duration FWHMs of 185 ± 10 and 270 ± 10 fs, respectively. Double-headed arrows in each plot indicate where the FWHM was measured. Fluorescent signal (Fl) is normalized to 1 at the center of the bead in (a), and at 1.5 ps in (b) and (c).

We characterized the temporal intensity distribution at the focal plane by incorporating a Michelson interferometer with a variable delay arm in the excitation path, and measuring the autocorrelation function of the focused excitation beam within a sample of 25-μM fluorescein. ^{36, 37} We calculated the pulse duration at the imaging plane from autocorrelations measured with the objective as a focusing lens and also when replacing the objective lens with a 1 in. diameter singlet lens with a 1 in. focal length. These two measurements [Figs. 2b and 2c] allow us to estimate the group delay dispersion (GDD) of the objective lens. We recorded the FWHM of a Gaussian fit, [TeX:] $\tau _{AC}^{{\rm FWHM}}$ ${\tau }_{AC}^{\mathrm{FWHM}}$ , to the moving average of our interferometric autocorrelation function over five measurements in each configuration. Assuming the original pulse shape is a Gaussian, the FWHM of the original pulse shape, [TeX:] $\tau _p^{{\rm FWHM}}$ ${\tau }_p^{\mathrm{FWHM}}$ , is related to the autocorrelation trace by: [TeX:] $\tau _p^{{\rm FWHM}} = 0.59\;\tau _p^{{\rm FWHM}}$ ${\tau }_p^{\mathrm{FWHM}}=0.59{\tau }_p^{\mathrm{FWHM}}$ , we estimate the pulse duration at the sample plane to be 185 ± 8 fs with the singlet lens and 270 ± 10 fs with the 20x/0.95 objective lens. Using the measured spectral bandwidth of 7.5-nm FWHM, and assuming a modest GDD of 250 fs² from the 1 in. singlet lens, we estimate that the GDD of our objective lens is approximately 4300 fs² at 760 nm. This GDD is significantly larger than reports of other similar NA objectives,³⁷ but not unexpected, given the long physical length of the lens (75-mm long).

2.3.

Monte Carlo Model

We used two independent Monte Carlo simulations to find the spatiotemporal distribution of the squared excitation intensity, I ²(r, z, t), and the spatial dependence of the collection efficiency, ϕ(r, z). These results were then combined using Eqs. 3, 8 to calculate R and determine image contrast, Q, as a function of imaging depth for a variety of sample and system parameters.

For the excitation simulation, the trajectories of photon packets were traced through a cylindrically symmetric three-dimensional grid with coordinates for depth z, radius r, and time t. All simulations used a Henyey–Greenstein scattering phase function with a scattering anisotropy, g, set to 0.85, which is the expected value from Mie theory for 760 nm light interacting with 1-μm diameter polystyrene spheres.³⁸ We used an absorption mean free path length of 5 mm in all simulations.

We modeled the focusing of the excitation light using the hyperbolic focusing method described by Tycho ³⁹ By skewing the initial trajectory of each photon off of a z-r plane by an angle depending on r, a focused Gaussian beam is exactly reconstructed in all three dimensions with photons that travel in straight paths. In a Gaussian beam focused at an imaging depth of z ₀, the lateral 1/e² beam radius as a function of depth, w(z), is defined by two parameters: the lateral 1/e² beam waist at the focal plane (w ₀) and the Rayleigh range (z _R):

Excitation photons were launched to a w ₀ equal to [TeX:] $\sqrt 2$ $\sqrt{2}$ times the 1/e² radius of the measured IPSF² at the focal plane. To find the z _R to match the profile of the excitation beam at the front aperture, we assume that the fill factor at the back aperture of the objective is linearly transferred to the front aperture of the objective—that is that w _BA/r _BA = w _FA/r _FA [Fig. 3a]. The resulting parameters describing the Gaussian distribution of the excitation photons are:

Fig. 3

(a) Schematic of parameters used in modeling signal and background fluorescence generation for a focused Gaussian beam in turbid media. (b) The ratio of signal fluorescence to background fluorescence in a homogeneously stained, R, is a function of the signal volume size, defined here as the volume of a sphere with radius r _s. The lines shown are exact analytical solutions and solid points are Monte Carlo results for focusing to a spot size equal to the measured ISPF² of our system (blue line), and for focusing to the conditions used in our model (red line) in the limiting case where the scattering and absorption of the media are zero.

We obtained the temporal impulse response of the system by setting initial photon times such that ballistic photons intercept the focal plane at exactly the same instant. The intensity distribution was determined by convolving the impulse response with the Gaussian pulse envelope measured in Sec. 2.2. Finally, the fluorescence generated per voxel was calculated by multiplying the integrated intensity squared by the volume of the voxel.

A second, time-independent, Monte Carlo simulation was used to model the effect of heterogeneous collection efficiency on maximum imaging depth. We implemented the approach outlined by Beaurepaire and Mertz, where the maximum radius and exiting angle of a photon escaping the tissue are coupled to conserve beam étendue.⁴⁰ We used an effective field of view radius of r _f = 1 mm.²⁷ We created a set of collection efficiency maps by launching photons with an isotropic initial angle distribution and a wavelength of 515 nm. The ratio of scattering length at 760 and 515 nm found to be 1.7 by Mie simulations for 1-μm polystyrene beads (see Sec. 2.5). Thus, for each excitation simulation set, we then used a scattering length of [TeX:] $l_s^{{\rm collection}} = l_s /1.7$ $l_s^{\mathrm{collection}}=l_s/1.7$ to create a matching collection efficiency map. We used an [TeX:] $l_a^{{\rm collection}}$ $l_a^{\mathrm{collection}}$ of 2 mm. To calculate R, the fluorescence grid created from the excitation simulation was integrated over time and every element was multiplied by an interpolated value from the spatial collection efficiency grid.

2.4.

Analytical Model

We applied the analytical model developed by Theer and Denk to calculate the fluorescence distributions for different optical system and sample parameters.⁴¹ In short, we consider the ballistic and scattered light distributions in turbid media independently and then combine the resulting intensities to calculate total fluorescence. The intensity from ballistic photons assumes an incident beam with Gaussian spatial and temporal distributions. The scattered light is modeled using a small angle approximation and statistical methods to calculate an effective temporal and transverse spatial beam width. We used an NA of 0.75, which produces a similar spot size to that measured from our IPSF². All other parameters used in the analytical model were the same as those used in the Monte Carlo model.

2.5.

Tissue Phantom Preparation

To measure the experimental contrast decay of two-photon imaging in a turbid media, we prepared 12 scattering agar phantoms with a range of optical properties similar to what is typically observed in epithelial tissues. Low melting point agarose (1.0%, Sigma) was prepared with 0.95-μm diameter polystyrene beads (Bangs Labs, PS03N) to vary scattering coefficient, 1-μm diameter fluorescent polystyrene beads (Invitrogen F-8823) to provide features for contrast measurements, and fluorescein (Fluka 46955) to control staining inhomogeneity. Polystyrene beads were ultrasonicated for 30 min to reduce aggregations before they were mixed into the agar phantoms. To increase the two-photon action cross-section of fluorescein, we mixed the agar solution with 2% pH 12 buffer. Polystyrene concentrations of 1.7×10¹⁰, 0.9×10¹⁰, and 0.6×10¹⁰ particles per mL provided scattering mean free paths of 40, 80, and 120 μm at 760 nm excitation, and 23, 46, and 69 μm at 515 nm emission, respectively. We determined the scattering coefficient with a Mie calculation using polystyrene sphere and agar gel refractive indices of 1.58 and 1.33, respectively.⁴² To verify the scattering lengths of our polystyrene solution matched Mie predictions, we measured attenuation of collimated 760-nm light through dilute solutions of polystyrene beads in a cuvette. Fluorescent beads were added at a constant concentration of 3.5×10⁸ particles per mL in all phantoms.

We added increasing amounts of fluorescein to get final concentrations of 0, 4, 10, and 25 μM for each scattering length tested, creating a total of 12 phantoms. This approach allows the formation of phantoms with low staining inhomogeneities without an increase in scattering coefficient that would result from adding high concentrations of fluorescent polystyrene beads. An additional benefit of using this method is that we can increase the background fluorescence immediately adjacent to the fluorescent beads. In contrast, previous studies that varied fluorescent bead concentration could control the fluorescence only in volumes many micrometers away from the signal volume.¹⁸ Thus our approach provides a more appropriate model for studying autofluorescence contrast decay, where the local staining inhomogeneity changes rapidly in short distances from the nucleus to cytoplasm to extracellular matrix. The disadvantage of this approach is that the out-of-focus background fluorescence cannot be measured directly because there is no region without some fluorescent signal—when the focal spot does not overlap with a fluorescent bead, the measured background still includes some contribution from the fluorescein in the focal volume.

We obtained lateral phantom images at 1 μm depth increments from the surface to the depth where the image contrast fell to approximately 0.1. The contrast for each of the 16 phantoms was analyzed using an automated bead-finding script written in MATLAB (Mathworks, Inc.). The MS for each bead was then calculated as the average of the maximum 4 voxels (150×150×1000 nm/voxel) within each region and the MB was calculated as the average pixel value of the nonconnected regions at the depth determined by the centroid of the bead.

2.6.

Biopsy Preparation

Human oral cavity biopsies were obtained from the University of Texas M. D. Anderson Cancer Center as part of an ongoing project for two-photon diagnosis of oral malignancies. The study was reviewed and approved by the internal review boards at M. D. Anderson Cancer Center and the University of Texas at Austin. Biopsies were excised from suspicious regions and contralateral normal tissue in the oral cavity and submitted for histopathology. We present two-photon images from the normal biopsy in this paper. The biopsy was approximately 3 mm in diameter and 5-mm thick, delivered in chilled culture media (Phenol Red-free DMEM High, Fisher Scientific), and imaged within 6 h of excision. The biopsy was stabilized on a Petri dish with low melting point agar and imaged at room temperature in a culture media bath.

3.1.

Monte Carlo Simulations

To validate our Monte Carlo model, we verified that for a transparent tissue, in the absence of scattering and absorption, the expected Gaussian profile is maintained from the sample surface to the full working distance of the lens. The value of R calculated from Monte Carlo simulations and the exact solution are shown in good agreement in Fig. 3b. Select Monte Carlo simulations were rerun with half the chosen grid spacings of r, z, and t (and corresponding eight times more grid elements for constant grid size), to verify the grid spacing was sufficiently small to capture the spatiotemporal changes in fluence. No significant changes were found in R calculations for smaller grid spacings, indicating the grid choice was sufficiently small to capture dynamic and spatial changes in intensity distribution. Finally, we verified that focal plane fluorescence generation scales inversely with pulse duration.

We used our Monte Carlo model to calculate R as a function of depth for each phantom. For these calculations, V _s was set to the volume of the 1 μm diameter fluorescent beads which provide the signal in the phantoms. Figure 4 presents the fluorescence distribution for imaging 400 μm deep in a homogeneously labeled sample with an l _s of 80 μm. In this case, the collected background fluorescence overwhelms the signal fluorescence by approximately 30 times. When the fluorescence is integrated circumferentially, we observe that the cumulative background fluorescence generation peaks slightly off-axis because the off-axis fluorescence is generated in larger volumes [Fig. 4a]. Finally, integration of this fluorescence radially shows that the total fluorescence generation at each transverse plane is relatively constant through the first two l _s, and monotonically decreasing at larger depths [Fig. 4b].

Fig. 4

Monte Carlo simulations are used to calculate the expected fluorescence distribution for an imaging depth of 400 μm in a sample with a 460-nm FWHM lateral spot size, τ_p = 270 fs, l _s = 80 μm, and g = 0.85. These parameters resulted in an R value of 1/30. (a) shows the lateral and axial distribution of the circumferentially integrated background fluorescence distribution. In this case, most of the background fluorescence comes from above the imaging plane and within 200 μm of the optical axis. (b) shows the contributions from ballistic (B ²), scattered (S ²), and combined 2(B+S) photons separately. The total fluorescence generation per transverse slice is relatively constant through the first three l _s deep. The solid and dotted lines represent simulations with excitation pulse durations of 270 and 135 fs, respectively. (a) is normalized so that the maximum out-of-focus fluorescence is one, (b) is normalized so the maximum fluorescence is one.

The contributions from ballistic (B) and scattered (S) fluence were separated by tagging scattered photons in the Monte Carlo simulations [Fig. 4b]. The fluorescence generation at each voxel will then be proportional to (S+B)². These results can provide some insight into how certain parameters affect the generation of background fluorescence. For instance, previous studies have shown that the R value can be increased by using shorter pulses. ^{18, 19} Looking at the individual components, as the pulse duration is decreased, it is entirely the scattered (S ²) and combined (2BS) terms that show a relative reduction in contribution of background fluorescence. Scattering can be thought of as a source of temporal dispersion from the sample on the scattered fluence. Furthermore, as the pulse duration is decreased, the dispersion from scattering reduces the fluorescence generated from scattered photons relative to the peak fluorescence generated within the focal volume, resulting in an increase in R for shorter pulse durations.

Though high collection efficiency is important for reaching large imaging depth with limited excitation power, we found that heterogeneous collection efficiency generally had a small effect on our calculations of R. In the l _s = 80-μm simulation ( [TeX:] $l_s^{{\rm emission}} = 46\;\mu {\rm m}$ $l_s^{\mathrm{emission}}=46\mu \mathrm{m}$ ), the collection efficiency only reduced R by 12% at z ₀ = 400 μm, and by 60% at z ₀ = 800 μm. The radial dependence of the collection efficiency had an especially weak effect on R, since the majority of the fluorescence generated in the excitation simulations was within 200 μm of the optical axis.

3.2.

Phantom Imaging

For all phantoms, we had sufficient power to reach contrast levels below one. Typically, excitation powers smaller than 1 mW were sufficient to image at the surface, while several hundreds of megawatts were required at the largest depths. For the phantom with l _s = 80 μm and χ = 300, for example, the excitation power was increased from 0.7 mW at the sample surface to 483 mW at the sample surface when imaging 510-μm deep. We identified an average of 190 beads per 77×77×100 μm³ field of view in each of the 12 phantoms using the automated bead-finding script. This number of beads corresponds to a bead concentration of 3.2×10⁸ beads per mL, slightly less than the expected value of 3.5×10⁸ beads per mL.

Figure 5 shows XZ images of the l _s = 80 μm phantom set for different staining inhomogeneities. Each XZ reconstruction is created from a maximum projection through 15 μm (45 pixels) of Y. The accompanied biopsy images are discussed in Sec. 3.5. We estimated the staining inhomogeneity of phantoms using the measured contrast at the surface of each phantom and the R values found from Monte Carlo simulations. The R value since it only changes slightly for shallow imaging depths in our Monte Carlo model—a decrease from 0.51 at the surface to approximately 0.46 at one l _s deep. Substituting an R value of 0.5 in Eq. 8, the staining inhomogeneity, χ, becomes approximately equal to three times the measured contrast at shallow depths, for large values of χ. Using this approximation, we estimate the average staining inhomogeneities to be χ = 300, 62, 25, and 10, by using the average contrast from beads identified within the first 20 to 40 μm of the three phantoms at each staining inhomogeneity.

Fig. 5

Comparison of XZ images of phantoms with constant scattering length of l _s = 80 μm for increasing staining inhomogeneity, χ, and a human tongue biopsy. Phantom cross sections are maximum projections through 15 μm of Y. Biopsy cross-sections shown are a maximum projection “max Biopsy” and a standard deviation projection “σ Biopsy” through 15 μm of Y. The standard deviation projection is normalized so that the maximum value is white.

The measured bead size remained constant with imaging depth for all phantoms. Figure 6 shows a representative case for the measured bead size versus depth for l _s = 80 μm, χ = 300 phantom. We determined the FWHMs of bead sizes by fitting a Gaussian function through the centroid of each bead found in the phantom. We found a slightly smaller bead sizes in the direction of our fast moving mirror (horizontal), indicating a slightly smaller resolution in that direction, possibly due to a slight ellipticity of the beam shape at the back aperture of the objective. The apparent axial resolution measured from the 1-μm beads was noticeably larger than that measured from the IPSF² with 100-nm beads. We attribute this increased measured axial resolution to sparse sampling (axial spacing between images was 1 μm) and long times between imaging the top and bottom of the beads (the time elapsed between the first image and the last image 3-μm apart was approximately 10 s). Nonetheless, the constant lateral and axial size of the measured beads does indicate minimal specimen-induced aberrations with increasing imaging depth. This result is in contrast to studies without index matching, where spherical aberrations are commonly observed. ^{43, 44}

Fig. 6

Measured sizes of 1 μm diameter fluorescent beads versus depth for l _s = 80 μm, χ = 300 phantom in the lateral and axial directions. The trend and error bars are calculated by the mean and standard deviations of sizes obtained by binning the beads at 50-μm depth increments. We observed no significant increase in bead size and, thus, in system resolution with increasing imaging depth in any of our phantoms.

3.3.

Fluorescence Decay

We measured the decay of the signal, background, and total fluorescence with increasing imaging depth in our 12 phantoms. By increasing the excitation power exponentially with imaging depth and normalizing the measured signal at each imaging depth by the expected quadratic increase with the excitation power, P(z = 0)²/P(z = z ₀)², we could accurately measure fluorescence decays far beyond the limited dynamic range of our detection system.

Figure 7 shows the measured and calculated fluorescence decays of the three χ = 62 phantoms. Looking at the χ = 62, l _s = 80 μm phantom, Fig. 7a shows the decays of the average value of the measured background at each imaging plane (MB), the measured signal for each fluorescent bead (MS), and the difference (MS-MB). Note that the slopes of the (MS-MB) and MB decays are approximately equal for imaging depths down to three mean free scattering lengths. This result suggests that the contrast, defined as the ratio of (MS-MB) to MB, is relatively constant for shallow depths. At larger imaging depths (z > 4l _s) the measured signal, which decays exponentially, is overcome by the background, which decays with [TeX:] $z_0^{ - 1}$ $z_0^{-1}$ .¹⁸ Comparing the fluorescence decay in phantoms with different scattering lengths, we found good agreement between the measured data and the Monte Carlo predictions [Figs. 7b and 7c ]. The decay of the fluorescent signal from the biopsy is discussed in Sec. 3.5.

Fig. 7

Plots of normalized fluorescence decays versus imaging depth for constant staining inhomogeneity (χ = 62). Points are data from phantom measurements and lines are the decays predicted by our Monte Carlo simulation with homogeneous (solid lines) and heterogeneous (dashed lines) collection efficiency. (a) The average fluorescence decay, 〈M〉, represents the average pixel value of a 512×512 pixel image recorded at each imaging depth. For comparison, the average fluorescence decay of the biopsy is also shown. (b) Examining relative decays for the MS, MB, and difference (MS-MB) versus depth for the χ = 62, l _s = 80 μm phantom, we found a maximum imaging depth of z _m = 390 μm. (c) The background-subtracted fluorescence decay exhibits exponential decay for the entire measured range. Monte Carlo simulations agree well with experiments for homogeneous collection efficiency (dashed lines) and heterogeneous collection efficiency (solid lines). Decays are normalized to one at 20-μm deep in (a) and (c). In (b), the background is normalized to one at 20-μm deep.

The effect of heterogeneous collection efficiency has little effect on the calculated decay curves, indicated by the overlap of the Monte Carlo results for shallow and moderate imaging depths. We also found that changing the staining inhomogeneity had little effect on the observed fluorescence decay rate, as summarized in Table 1. As expected, the staining inhomogeneity had the effect of raising or lowering the initial MB value relative to the MS values.

Table 1

Summary of fluorescence signal decay length constants from exponential fits to measured phantom data (from z 0 = 0 to 3l s) and predicted in Monte Carlo data, with and without including the effect of heterogeneous collection efficiency.

Table 1 summarizes the measured length constants, [TeX:] $l_s^m$ $l_s^m$ , obtained from the exponential decay rates of the background-subtracted fluorescence in each phantom, and the corresponding decay constants predicted by Monte Carlo models. In all phantoms, the actual scattering length could be estimated to within 10% accuracy by using the approximation: [TeX:] $l_s = 2.3l_s^m$ $l_s=2.3l_s^m$ . This result is reasonable considering a previous report, which found a factor of 2.5 relation between [TeX:] $l_s^m$ $l_s^m$ and l _s using an NA of 1.2.¹⁵

3.4.

Contrast Decay

The measured contrast decay from tissue phantoms showed similar trends to those predicted by Monte Carlo and analytical models (Fig. 8). The Monte Carlo model matches the analytical model relatively well for the higher χ values tested. For the χ = 10 phantom, the measurements demonstrate shorter z _m than predicted, likely due to the high concentration of fluorescein in these phantoms beginning to absorb significant amounts of emission (the “inner filter” effect). If the absorption mean free path length is known at the emission wavelength, this effect could be accounted for using a shorter mean free absorption length in our collection efficiency calculations.

Fig. 8

Contrast decays of phantoms with (a) constant staining inhomogeneity and (b) constant scattering length. Monte Carlo simulations (solid lines) agree well with the analytical model (dashed lines) and the phantom contrast measurements (solid dots). Both models slightly overestimate the maximum imaging depth, increasingly at lower staining inhomogeneities.

We observed that for the χ = 300 phantoms, the contrast was higher than the expected values at shallow depths [Fig. 8b]. We attribute this large contrast to an effective increase in χ at shallow depths. Because no fluorescein was added in this set of phantoms, the staining inhomogeneity was determined by the influence of background fluorescent beads. For the shallowest beads, there are no background fluorescent beads directly above them, effectively increasing their apparent staining inhomogeneity. For deeper beads, the fluorescent beads that are present above the imaging plane produce background fluorescence.

3.5.

Human Biopsy Imaging

Figure 9 presents a summary of the 200 autofluorescence images collected at 2-μm depth increments from a biopsy of healthy human tongue tissue. The water-dipping objective allowed us to visualize the natural surface roughness of the tissue, which plays a role in collection efficiency⁴⁰ [Fig. 9a]. We defined the surface (z = 0 μm) as the depth at which half of the field of view had signal. We increased the excitation power delivered to the tissue surface gradually from 3 mW at the surface to 30 mW at 170-μm deep. At larger imaging depths, we could maintain constant signal detected while increasing the excitation power less rapidly. We used a maximum of 160 mW of excitation power delivered to the sample surface when imaging 380-μm deep.

Fig. 9

(a) Three-dimensional rendering of a sequence of 200 lateral 2PAF images acquired from healthy human tongue biopsy. (b) Selection of lateral images from an imaging depth of 40 to 360 μm. Field of view in all lateral images is 170 μm. (c) Normalized signal profile from manually identified cells at imaging depths of 40, 120, 240, and 360 μm. Profiles are taken from lateral lines indicated in (b).

Lateral images show subcellular resolution and cellular morphology with bright cytoplasm to dark nucleus contrast [Fig. 9b]. Contrast is presumably due to high concentrations of NADH, NAD, and FAD in the cytoplasm, which have been shown to be the dominant endogenous fluorophores in confocal autofluorescence imaging of cervical and oral epithelial tissues, ^{45, 46} and are efficiently excited in 2PAM at 760-nm excitation.⁴ We found especially high signal levels from the stratum corneum, which, based on our Monte Carlo model, would strongly contribute to the out-of-focus fluorescence found at large imaging depths.

A maximum XZ projection through 15 μm of Y shows that at large depths, bright pixels in the background substantially degrade the contrast beyond 300 μm [Fig. 5]. A normalized standard deviation projection through 15 μm of Y gives values similar to the measured contrast—it shows bright regions where there are large changes in pixel values. At shallow depths, this projection shows large values where there is variation between the bright cytoplasm signal and dark nuclei. At large depths, the increase in low-spatial-frequency, out-of-focus background fluorescence reduces the lateral variation in signal, and the corresponding standard deviation along Y decreases.

We examined the total fluorescence decay with imaging depth to estimate the scattering coefficient of the biopsy. Though fluorescence decay accurately predicted the scattering coefficient in phantoms, extending this method to the biopsy can be complicated by several effects. 1. The scattering length and fluorophore concentrations are likely to change with depth. Previous studies have observed brighter signal coming from the stratum corneum and basal membrane than at the intermediate epithelial layers. ^{46, 47} 2. As the excised tissue dies, the autofluorescence signal levels gradually decay with time. We observed that the superficial regions of the tissue begin losing fluorescence signal more rapidly than the interior regions, creating a time-dependent change in relative fluorescence density. 3. Specimen induced aberrations typically become more severe with increasing imaging depths,⁴⁸ and would result in a higher apparent rate of fluorescence decay. Though we verified aberrations played a minimal role in degradation of signal in the phantom, without a constant-size point source in our biopsy, we were unable to monitor our excitation spot size with depth in the biopsy. Nevertheless, we can obtain an estimate of the biopsy scattering coefficient by comparing our measured biopsy fluorescence decay to phantom experiments [Fig. 7a]. We fit an exponential curve to the first 200 μm of fluorescence decay obtained from the biopsy and found a length constant of 39 μm, corresponding to a scattering length of 89 μm, assuming [TeX:] $l_s = 2.3l_s^m$ $l_s=2.3l_s^m$ . At larger depths, the fluorescence decayed less rapidly, showing a similar slope to the l _s = 120-μm phantom.

Figure 9c presents typical contrast levels obtained at different imaging depths. We examined the contrast in line profiles drawn through manually identified epithelial cells marked by the lines drawn in the lateral images displayed in Fig. 9b. At 40-μm deep, the brightest signal from the cytoplasm is approximately 20 to 40 times larger than the minimum signal from within the nucleus, while the average signal from the cytoplasm is approximately 5 to 10 times the average signal from the nucleus [Fig. 9c]. Using averaged signal and background as contrast, and employing the same method used with the phantoms, we estimate the χ of the biopsy to be 15 to 30. The contrast falls to 1 in our biopsy around 320 μm below the tissue surface, approximately 3 to 4 mean free scattering lengths deep.

3.6.

Fluorescence Saturation and Photobleaching

The analysis presented in this paper has assumed that the intensity ranges used in our experiments are low enough to assume negligible fluorescence saturation and photobleaching—that is, that the fluorescence generation is quadratically dependent on the excitation intensity and is time-independent. However, at high intensities, these assumptions no longer hold. We performed a series of tests to verify that we do not experience any saturation and photobleaching in our experiments.

The average powers used for surface imaging in the phantoms and biopsy were 1 to 3 mW, respectively. Using our measured spot size and pulse duration, these powers correspond to peak intensities of 5 to 15 GW/cm² generated at the focal spot. These intensities and the estimated ∼320 overlapping pulses per spot used in our study, are well below imaging conditions found in the literature to have no affect cell viability, as summarized by Hoy ⁴⁹ Furthermore, we found in phantom experiments and Monte Carlo simulations that the rate of power increase necessary for maintaining constant fluorescence detected is close to that required to deliver constant fluence to the focal point. Though when imaging deep inside scattering media, average powers greater than 100 mW are delivered to the sample surface, it is unlikely that the fluence delivered to the imaging plane at large depths is significantly greater than that delivered to the imaging plane when imaging the sample surface.

Fluorescence saturation can broaden the point spread function, as fluorophores in the highest intensity regions no longer generate a greater emission signal than fluorophores further away. The resulting decrease in resolution has been previously demonstrated, but is expected to have a significant effect only under a combination of high intensities and especially large two-photon action cross-section⁵⁰ (>1000 GM cross-section at 1 mW of excitation power delivered to the imaging plane, where 1 GM unit equals 10⁻⁵⁰ cm⁴ s). Our observation of constant measured bead size with increasing imaging depth indicates a negligible influence of fluorescence saturation in our phantom experiments [Fig. 6]. To test for fluorescence saturation in our biopsy experiments, we verified that our detected signal scales quadratically with excitation at select imaging depths. We measured power dependencies of 1.94, 2.00, 1.97, and 1.99 at increasing imaging depths of z ₀ = 290, 300, 326, and 356 μm, respectively. In phantoms, we also observed that the signal from the identified beads scaled to the power of 1.98±0.05 at each 100 μm imaging depth increment.

We also tested for the presence of photobleaching in 2PAM of the biopsy by measuring the ratio of the average signal of the first and last image taken at every image plane. Each imaging plane was raster scanned 8 times at 3 frames per second. We found an average ratio of 1.00 ± 0.02 for the biopsy, and no trend for photobleaching at larger imaging depths. In conclusion, we expect that fluorescence saturation and photobleaching did not appreciably influence the results presented herein.

3.7.

Maximum Imaging Depth

Figure 10a summarizes the maximum imaging depth found in phantom experiments, Monte Carlo simulations, and the analytical model for samples with different scattering lengths and staining inhomogeneities. We found reasonable agreement between the three approaches, with experiments matching the models slightly better for longer scattering lengths and higher staining inhomogeneities. This difference is likely due to the increasing influence of small heterogeneities in our shorter scattering length phantoms, such as spatial variations in fluorescein and polystyrene bead concentrations. The Monte Carlo model showed a slightly weaker dependence of maximum imaging depth on scattering length, partially due to the effect of inhomogeneous collection efficiency.

Fig. 10

(a) The maximum imaging depth determined by the depth at which Q = 1. (b) Expressing the maximum imaging depth in terms of scattering mean free paths, we observed a linear dependence of maximum imaging depth on log(l _s). Data are plotted from phantom measurements, as well as analytical and Monte Carlo models, with homogeneous collection, and real heterogeneous collection.

The approximate position of our biopsy is also indicated, based on extraction of optical properties from biopsy images. The maximum imaging depth achieved in the biopsy was 320 μm, which was approximately 20% less than that expected from phantoms with similar l _s and χ. This difference is likely due to greater specimen-induced aberrations in biopsy imaging.

Normalizing the maximum imaging depth by [TeX:] $z_{m}^ n=z_m/l_s$ $z_m^n=z_m/l_s$ , we find a logarithmic dependence on l _s [Fig. 10b]. Though this is a reasonable approximation for small fractional changes in l _s, it leads to an overestimate of maximum imaging depth by 1 scattering length when extrapolating maximum imaging depths found in samples with l _s = 200 μm to samples with an l _s of 40 μm. The origin of the dependence of z _m on l _s is that as the distances are scaled down by l _s, the photons responsible for the background fluorescence are confined to smaller volumes while the photons generating focal volume signal pass through a constant size signal volume.