2.1.

Experimental Setup

A home-built, inverted laser scanning microscope was used to perform TPM and THG microscopies [Fig. 1(a)]. The excitation laser was a femtosecond Er-doped fiber laser (Discovery, Raydiance Inc.) with 3-W average power at 1552 nm ($1.5\mu \mathrm{J}$ pulse energy) and 1 W when frequency doubled to 776 nm ($0.5\mu \mathrm{J}$ pulse energy) at 2 MHz with 600-fs pulse width. The laser beam was focused with a 0.75 NA, $20\times$ air objective (Nikon Plan Apo). We measured the spot size by imaging 100-nm-diameter polystyrene beads (Invitrogen, F8803) embedded in agar gel at depths varying from 50 to $500\mu \mathrm{m}$. The agar-bead phantom consisted of $10\mu \mathrm{l}$ of the stock bead solution (2% solids) dispersed in 2 ml of 4% agar in water, which was then allowed to gel. This protocol resulted in a bead concentration of $2\times {10}^{10}\text{beads}/\mathrm{ml}$, which was estimated to yield an $\sim 35\mu \mathrm{m}$ scattering length at 776 nm excitation wavelength, consistent with the effective scattering lengths expected from the lamina propria (LP) of our samples. The average lateral full-width half-maximum (FWHM) of the two-photon point spread function (PSF) was $0.58\pm 0.07\mu \mathrm{m}$ and three-photon FWHM of the PSF was $0.95\pm 0.11\mu \mathrm{m}$. These resolutions correspond to a $1.39\pm 0.17$ and $2.8\pm 0.26\mu \mathrm{m}$ $1/e^2$ diameter of the intensity distribution at 776 and 1552 nm, respectively.

Fig. 1

Schematic of the inverted microscope system for nonlinear imaging and ablation. (a) Femtosecond laser pulses were obtained from a compact fiber laser system at the fundamental wavelength (1552 nm) and passed through an energy attenuator consisting of a half-wave plate ($\lambda /2$) and polarizing cube beam splitter (PCBS). Frequency doubled pulses (776 nm) were achieved using a beta barium borate crystal. Both beams were scanned by a pair of galvanometric scanning mirrors (SM), which was then imaged on the back aperture of a 0.75 NA, $20\times$ objective by a pair of a scan lens (SL) and a tube lens (TL). The samples were placed on a three-axis motorized stage (P280, Nanocube, PI) for nonlinear imaging. Emitted light was collected by a dichroic mirror (DM), collection optics (CO), laser blocking filters (Edmund Optics, 84-650 for THG imaging, or Semrock, FF01-750/SP-25 for TPM imaging), and either a THG filter (Chroma, HQ515/15M) or a TPM filter (Schott, BG38) into the photomultiplier tube (PMT) (H10770 PA-40, Hamamatsu). (b) The slopes of signal intensities versus imaging power verified the order of the nonlinear processes for both TPM and THG imaging on a log-log scale plot.

To satisfy the Nyquist criterion for both TPM and THG, a pixel size of $0.3\mu \mathrm{m}$ was used corresponding to half of the resolution of TPM ($0.6\mu \mathrm{m}$). The laser beam was scanned in the $x$ and $y$ axes using a pair of galvanometric mirrors (Cambridge Technologies, Inc.) and sweeping the focal spot over a $150\times 150{\mu \mathrm{m}}^2$ field of view for 3.3 s for both modalities to obtain an image of $512\times 512$ pixels. The resulting number of pixels corresponded to 2 and 3 pulses/pixel/raw image for TPM and THG imaging, respectively. We collected and averaged 10 images at each plane. A dichroic mirror (Semrock, FF735-Di01) was used to transmit the excitation light to the sample and reflect the emitted signal from the sample to the collection optics. A photomultiplier tube (PMT) (Hamamatsu, H10770 PA-40) collected the signal via collection optics, laser blocking filters (Edmund Optics, 84-650 for THG imaging, or Semrock, FF01-750/SP-25 for TPM imaging), and emission filters (Chroma, HQ515/15M for THG imaging and Schott, BG38 for TPM imaging). The computer software (MPScan) reconstructed these signals into $512\times 512$ pixel images at a frame rate of 3.05 Hz.³⁶

The order of nonlinear processes was characterized by measuring the detected signal on the PMT while increasing laser power delivered to the sample. As shown in Fig. 1(b), the slope of the PMT signal versus imaging power curve yielded $3.02\pm 0.12$ and $2.05\pm 0.08$ on a log-log scale for THG and TPM, respectively. As expected, these results confirmed that THG and TPM imaging were second- and third-order nonlinear processes.

2.2.

Ex Vivo Tissue Preparation

The LP of porcine vocal folds has been previously studied and was found to be most similar to the human LP in terms of their elastin and collagen distribution.^37,38 Our previous study also showed that superior porcine vocal folds have a fiber structure similar to human vocal folds as compared to inferior porcine vocal folds at the depths we are interested in.¹⁰ Thus, it was decided to use freshly excised superior porcine vocal folds in the present experiments. Fresh porcine airway specimens were acquired from a local slaughterhouse and the larynx was isolated in a room temperature saline bath within 2 h of sacrificing the animal. After excision, the vocal folds were placed in saline solution and covered with a glass cover slip to flatten their upper surface. In a clinical application, the glass cover slip would be akin to the window of a microsurgery probe in contact with the sample, thus maintaining a constant depth of ablation. Short-lived autofluorescence of fresh tissue samples and fluorescence beads helped identifying the surface of the sample during TPM imaging. For THG imaging, high contrast at the air and glass interface allowed to identify the tissue surface easily.

3.1.

Extinction Properties of Tissue

To determine extinction properties of fresh porcine vocal folds at both excitation wavelengths (776 and 1552 nm), a novel ablation-based method was used that we developed in recent work.¹⁰ In this method, ablations are performed at three different depths below the epithelium in the SLP (72, 90, and $108\mu \mathrm{m}$) and the minimum pulse energy required to initiate ablation ($E_{\mathrm{th},\text{surface}}$) is found at each depth. For these ablation measurements, the laser is used at 303-kHz repetition rate. Beer’s law can be applied for each ablation depth, $z_{\mathrm{ab}}$, which was chosen to be below the measured epithelium thickness $z_{\mathrm{ep}}$ as

For 776-nm excitation wavelength, the tissue extinction length is equal to the scattering length because of the negligible tissue absorption. However, the absorption is considerable at 1552 nm. To identify the scattering and absorption components of the extinction length at 1552 nm, the scattering length is first estimated using the following empirical correlation for fibrous soft tissues:⁴¹

The resulting scattering length was $108\pm 8\mu \mathrm{m}$ at 1552-nm wavelength according to Eq. (4) and the absorption length could then be deduced as $135\pm 15\mu \mathrm{m}$. Several studies have found that the absorption length in pure water is $\sim 200$ to $250\mu \mathrm{m}$ at 1552-nm wavelength.^42,43 To verify these numbers, the absorption length of distilled water was also measured by focusing 1552-nm laser light in different depths and measuring the minimum pulse energy that creates bubble formation at the focal plane. By applying the method described above, this study determined the absorption length at 1552-nm wavelength to be $200\pm 50\mu \mathrm{m}$ in agreement with the literature. Thus, we believe that these tissue experimental results represent relevant numbers for vocal folds, and the small difference in the values of absorption lengths between water and vocal folds can be attributed to the vocal fold constituents such as collagen, elastin, and hyaluronic acid in addition to water. These values were used in the Monte Carlo and heat transfer simulations. The summary of all extinction lengths, which were generally in agreement with the values reported in the literature^41,44 for both excitation wavelengths, are tabulated in Table 1.

Table 1

Summary of measured tissue optical properties at 776- and 1552-nm excitation wavelengths. Extinction lengths for epithelium (ℓext,ep) and lamina propria (ℓext,LP) at the maximum imaging depths for two-photon and third-harmonic generation imaging at the corresponding wavelengths. The effective scattering (ℓsca) and absorption lengths (ℓsca,eff) at the maximum imaging depths (zm) are also presented. Standard deviations are included for each parameter as derived from the variations in the measured epithelial thickness.

3.2.

Maximum Imaging Depth

An excised fresh superior porcine vocal fold was imaged using both TPM and THG imaging to compare their maximum imaging depths at 776- and 1552-nm excitation wavelengths, respectively (Fig. 2). The tissue surface was identified by the presence of externally applied fluorescent beads in TPM and by the presence of large signal contrast at the interface between the tissue surface and the glass cover in THG. The depth at which collagen fibers first appeared varied between the locations of imaging, ranging from 40 to $70\mu \mathrm{m}$ below the surface. As the imaging depth was increased, TPM and THG signals were detected originating from different tissue components, which helped to differentiate epithelium and SLP. For instance, TPM at 776-nm excitation wavelength could detect nicotinamide adenine dinucleotide phosphate and flavins in the epithelium,³ and collagen and elastin fibers in the SLP^45,46 [Fig. 2(b)]. On the other hand, THG allowed to image epithelial cells and ECM fibers in superior porcine vocal folds. The THG signal originated from the index of refraction changes, geometric shape of scatterers, and third-order susceptibility differences inside the tissue [Fig. 2(a)]. Comparison between the nonlinear images of a superior porcine vocal fold clearly showed that THG at a longer wavelength (1552 nm) qualitatively improved the maximum imaging depth substantially. At the maximal imaging depths, the applied maximum imaging power was 10 and 55 mW for TPM and THG, respectively.

Fig. 2

Representative nonlinear optical images of a fresh superior porcine vocal fold. (a) Images of THG showing clear collagen features beyond $400\mu \mathrm{m}$ using imaging powers in the range of 50 mW. (b) Images of combined TPM and second-harmonic generation indicating maximum imaging depths limited to below $150\mu \mathrm{m}$. The scale bar represents $50\mu \mathrm{m}$.

To quantitatively analyze the maximum imaging depths, we analyzed how the SBR varied with imaging depth (Fig. 3). The SBR was calculated by selecting a line of pixels across resolvable features at different locations in each image using Image J software. These resolvable features were either epithelial cells in the epithelium or extracellular matrix fibers in the SLP. Then, the ratio was determined of maximum signal intensity to the average of 10 pixels on the right and left tail of the $1/e^2$ of the maximum signal intensity. The SBR for each imaging modality is presented with respect to the imaging depth on a semilog plot in Fig. 3(a). The maximum imaging depth, $z_m$, was defined as the depth where the SBR fell to 1. In addition, SBR was plotted against imaging depth normalized with tissue extinction length, ${{\ell }}_{\mathrm{ext}}$ [Fig. 3(b)].

Fig. 3

Signal-to-background ratio (SBR) values for both nonlinear imaging modalities: (a) SBR with respect to imaging depth, $z_m$, and (b) SBR with respect to normalized imaging depth, $z_m/{\ell }_{\mathrm{ext}}$.

Figure 3(a) quantitatively shows that THG imaging provided three times deeper imaging, increasing the imaging depth from $140\mu \mathrm{m}$ (TPM at 776 nm) to $420\mu \mathrm{m}$ (THG at 1552 nm). Figure 3(b) also shows that SBR became one at normalized imaging depths of approximately four and seven extinction lengths for the TPM and THG imaging modalities, respectively. Such improvement of nondimensional depths ($\sim 1.7$ times) was a direct result of the reduced probability for out-of-focus generation in three-photon process as well as the reduced absorption of the signal generated at longer wavelengths. Overall, both longer excitation and emission wavelengths and three-photon processes had improved the maximum imaging depth for THG microscopy at 1552-nm excitation wavelength. This maximum imaging depth improvement reveals a potential application in imaging scarred human vocal folds, which have thicker epithelium compared to porcine vocal folds.

3.3.

Analytical and Monte Carlo Modeling of Maximum Imaging Depths

To model the maximum imaging depths, analytical and Monte Carlo methods used previously by our group were used for varying extinction lengths and inhomogeneities ($\chi$)⁵. Assuming the concentration of fluorophores within focal volume is constant, $C_s$, and that the out-of-focus fluorophore concentration is diffuse enough to be approximated by the average fluorophore concentration, $C_b$, $\chi$ was defined as the ratio of $C_s$ to $C_b$. A range of $\chi$ values (10, 25, 62, and 300) has been chosen in the simulations to broadly represent various tissue types. The experimental results, represented within the circles with dashed lines in Fig. 4, are based on the experimental extinction lengths and normalized imaging depths. Interestingly, the $\chi$ value of 32, which was extracted from the THG images in Fig. 2, is in agreement with the $\chi$ values of 25 to 62 present within the circle that represents the THG experimental results [Fig. 4(b)].

Fig. 4

Analytical and numerical modeling of maximum imaging depths: (a) for TPM at 776-nm excitation and 515-nm emission wavelengths and (b) for THG at 1550-nm excitation and 515-nm emission wavelengths. The maximum imaging depths were calculated as a function of the effective extinction length (${\ell }_{\text{ext}}$) and for different tissue inhomogeneities ($\chi$). The extinction length varied for different scattering lengths representing different tissue types and assuming a negligible absorption length at 776 nm and ${\ell }_{\text{abs}}=135\mu \mathrm{m}$ at $1552\mathrm{nm}$. THG results were presented for extinction lengths smaller than $120\mu \mathrm{m}$, corresponding to a maximum scattering length of $600\mu \mathrm{m}$ (expected values for brain tissue at 1552 nm). Cross marks represent Monte Carlo simulation results and circles with dashed lines represent the approximate range of our experimental results.

This study adapted the analytical method developed by Theer and Denk⁶ to calculate the ballistic and scattered light distributions in turbid media, which were then combined to obtain the total fluorescence distribution with depth. The intensity for ballistic photons was assumed to have Gaussian spatial and temporal distributions. Monte Carlo simulations were used to simulate photon transport in the tissue with an anisotropy factor of $g=0.85$ by separately modeling the excitation and emission phenomena. Simulations calculated the spatiotemporal distribution of the intensities ($I^2$ and $I^3$) and the spatial dependence of collection efficiency. The image contrast could then be calculated as the ratio of the signal in the focal volume to the background. The maximum imaging depth was defined as the depth at which the image contrast fell to 1. In the case of THG, the emitted photons were directed forward, unlike the two-photon autofluorescence process where emission was assumed to be isotropic.^47,48

Figure 4 presents the results from Monte Carlo simulations and analytical calculations for different extinction lengths. For TPM imaging, normalized imaging depth had a logarithmic dependence on extinction length, as expected according to the Beer-Lambert’s law [Fig. 4(a)]. Monte Carlo simulations correlated well with the analytical calculations. The normalized imaging depth was estimated to be between four and five extinction lengths for a $\chi$ range of 25 to 62 for the range of scattering lengths expected in tissues (which is equal to extinction since absorption is minimal at 776 nm). These calculations were consistent with the experimental results shown in dashed circular box in Fig. 4(a).

Similar to the TPM case, normalized maximum imaging depth in THG microscopy also exhibited a logarithmic dependence on extinction length and increased with the increase in $\chi$ for all extinction lengths [Fig. 4(b)]. The measured extinction length in the experiments was $60\mu \mathrm{m}$ and the corresponding maximum imaging depth was estimated to be in the range of 7.0 extinction lengths. This result closely matched the simulated results for a $\chi$ range of 25 to 62. Another important conclusion drawn from Fig. 4(b) was that THG microscopy improved maximum imaging depth approximately two times compared to TPM microscopy at the same extinction lengths for all inhomogeneity values. This conclusion could be attributed to the effect of three-photon processes on reducing background signal and thus increasing SBR value.

3.4.

Tissue Heating Properties During Third-Harmonic Generation Imaging

To avoid thermal damage during THG imaging, the heating properties of tissue were investigated by measuring tissue surface temperature using an infrared thermal camera (FLIR, A325SC). Heat transfer simulations also were performed to predict the temperature distribution below the surface. To keep the field of view clear for the thermal camera’s view, the laser light was focused on tissue surface using a long working distance objective (Nikon, Plan Apo, $10\times$, 0.3 NA) for surface temperature measurements.

From these surface temperature measurements, the study first characterized thermal relaxation time of the tissue, $\tau$, defined as the time that the temperature dropped to $1/e$ of its maximum value after the laser was shut down. Knowledge of this parameter enabled to perform THG imaging without causing substantial heating of tissue by blocking the laser between different imaging planes for the duration of the characteristic cooling times. The thermal relaxation time of the bulk tissue was measured for different laser excitation conditions where the tissue was juxtaposed against a cover glass, and the laser was scanned for 4 s while being focused at the maximum imaging depth of $420\mu \mathrm{m}$ within the tissue. The measured temperature distribution for a single imaging session [Fig. 5(a)] showed a maximum temperature increase of 1 to 6°C as the imaging power was increased from 120 to 480 mW. The thermal relaxation times ($\tau$), tabulated in Table 2, varied only slightly around 15 s, which were in the range of the values cited in the literature for various tissue types.^44,49 The average value was also in agreement with an estimated value of 14 s that could be obtained using the simplified expression for relaxation time

Fig. 5

Tissue surface temperature measurement using an IR camera. (a) Measurement of tissue surface temperature to calculate tissue relaxation time for three imaging powers applied for the duration of 4 s at the maximum imaging depth of $420\mu \mathrm{m}$. (b) Comparison of measured surface temperature evolution with our heat transfer simulation results and the laser power used at a given depth. Experimental results correlated well with the experiments for a thermal relaxation constant of 15 s. The laser was focused inside tissue using a long working objective lens to clear the imaging path for the IR camera.

Table 2

Tissue thermal relaxation times (τ) for exposure to different average powers (Pave).

Since the maximum temperature was expected to be close to the focal volume below the surface where the temperature could not be measured, heat transfer simulations were performed to estimate the depth-resolved temperature distribution. The numerical analysis first calculated the intensity distribution inside the tissue by using the transient radiative transport equation (RTE).^50,51 It then calculated the temperature distribution by coupling Pennes energy equation with a non-Fourier damped heat conduction equation. ^50,51 It has been previously shown that the non-Fourier hyperbolic heat conduction model was a more accurate model since it takes into account the relaxation time of the tissues. In addition, the Fourier parabolic heat conduction model is found to deviate significantly from experimental measurements while studying heat transfer phenomenon for time scales shorter than thermal relation time of the medium.⁴⁴

To validate the numerical results, we measured the evolution of surface temperatures and imaging power over the duration of a depth-resolved THG imaging [Fig. 5(b)]. To avoid heating the tissue, the laser was blocked in between different imaging planes (separated by $12\text{-}\mu \mathrm{m}$ steps) for two thermal relaxation time periods (30 s). The measured surface temperature distribution correlated best with the numerical temperature distribution assuming a 15-s thermal relaxation time in agreement with the measured data. The corresponding results indicated that the maximum surface temperature did not exceed 32°C, representing only a 5.5°C increase from the initial room temperature. These results were obtained with a long working distance objective for an average excitation power ranging from 1 mW at the surface up to 100 mW at the maximum imaging depth of $420\mu \mathrm{m}$. This result verified both the experimental and numerical methods to determine depth-resolved tissue temperature distribution for THG imaging.

After this verification, numerical analysis was applied for a typical 3-D THG imaging where the high NA objective (0.75 NA, $20\times$, Nikon) was used and the excitation power was varied with depth as indicated in Fig. 2. The estimated surface temperature and imaging power over the duration of this experiment are shown in Fig. 6(a). Figure 6(b) represents estimated depth-resolved temperature distribution. A maximum temperature increase of $\sim 2°\mathrm{C}$ occurred at the maximum imaging depth of $420\mu \mathrm{m}$, where imaging power was maximum (55 mW). This temperature increase was not deemed detrimental to the tissue constituents.^52,53 This temperature model was based on linear absorption. To elucidate the effect of nonlinear absorption on the transient temperature evolution at the focus, the models of Vogel and his colleagues were used.⁵⁴ This nonlinear absorption analysis resulted in $\sim 0.3°\mathrm{C}$ temperature increase for both a single pulse and series of pulses. Thus, nonlinear absorption was safely neglected in the presented thermal model. This result showed how accurately one can control temperature increase by measuring the thermal relaxation time and manipulating the laser shuttering time in THG imaging.

Fig. 6

Dynamics of tissue temperature during typical THG imaging conditions. (a) The evolution of the surface temperature and imaging power as we image deeper into tissue. (b) Axial distribution of tissue temperature as modeled using a 15-s time constant. The maximum temperature increase was observed at the central line of the focal plane at the maximum imaging depth ($420\mu \mathrm{m}$) with maximum imaging power (55 mW) and was $\sim 2°\mathrm{C}$. (c) Tissue focal plan temperature and experimental duration versus interval time between consecutive THG imaging with a 0.75-NA objective lens.

Finally, the maximum temperature was estimated inside the tissue at the focal plane for different designated periods of time between consecutive imaging planes separated by $12\text{-}\mu \mathrm{m}$ steps [Fig. 6(c)]. An interval time of 30 s (twice the thermal relaxation time) indeed provides an optimized condition from the point of view of both experimental duration and temperature increase in the focal plane. For these conditions, the experimental duration is $<20\min$ and maximum temperature increase is $<2°\mathrm{C}$.