Imaging deeper than the transport mean free path with multiphoton microscopy.

Multiphoton fluorescence microscopy enables deep in vivo imaging by using long excitation wavelengths to increase the penetration depth of ballistic photons and nonlinear excitation to suppress the out-of-focus fluorescence. However, the imaging depth of multiphoton microscopy is limited by tissue scattering and absorption. This fundamental depth limit for two-photon microscopy has been studied theoretically and experimentally. Long wavelength three-photon fluorescence microscopy was developed to image beyond the depth limit of two-photon microscopy and has achieved unprecedented in vivo imaging depth. Here we extend the theoretical framework for characterizing the depth limit of two-photon microscopy to three-photon microscopy. We further verify the theoretical predictions with experimental results from tissue phantoms. We demonstrate experimentally that high spatial resolution diffraction-limited imaging at a depth of 10 scattering mean free paths, which is nearly twice the transport mean free path, is possible with multiphoton microscopy. Our results indicate that the depth limit of three-photon microscopy is significantly beyond what has been achieved in biological tissues so far, and further technological development is required to reach the full potential of three-photon microscopy.


Introduction
Optical microscopy plays a key role in advancing our understanding of biological systems. The high spatiotemporal resolution of noninvasive optical imaging techniques make them unique in their ability to answer specific biological questions [1][2][3][4]. Multiphoton microscopy (MPM) enables deep imaging within turbid media with high spatial resolution [5][6][7]. The nonlinear nature of multiphoton excitation reduces out-of-focus (background) fluorescence generation. Additionally, by using longer wavelengths to excite widely available fluorophores, the excitation photons are scattered less and can form a sharp focus deeper inside turbid media [5,8]. Over the last few decades, two-photon microscopy (2PM) has become the gold standard of deep in vivo imaging, enabling access to depths beyond the superficial layers [9][10][11][12]. While 2PM can penetrate deeper than one-photon confocal fluorescence microscopy when imaging the same fluorophores, it is eventually limited by contrast degradation due to background fluorescence generation [8]. As powerful ultrafast lasers have become more accessible, the depth limit of 2PM has been reached and quantitatively characterized [8,12]. For example, both theoretical and experimental investigations in the past have shown that the fundamental depth limit for 2PM is approximately 0.8 mm penetration depth when imaging mouse brain vasculature in vivo at the excitation wavelength of 920 nm [8,13].
To overcome the depth limit of 2PM, long wavelength three-photon fluorescence microscopy (3PM) was developed. 3PM takes advantage of the optimal excitation wavelength windows for deep tissue imaging at 1300 nm and 1700 nm and further suppresses background fluorescence because of higher order nonlinear excitation [5,14,15]. 3PM has since achieved unprecedented imaging depths; examples include 3PM of neuronal activity with single-cell resolution in the mouse hippocampus [6] or through the intact skull of adult mice and zebrafish [7,16], imaging mouse brain vasculature and blood flow at >2 mm depth [17], imaging the physiological response of the subplate of the mouse visual cortex [12,18,19], two-photon (2P) and three-photon (3P) volumetric imaging for simultaneous monitoring of both the cortical layers and the subcortical regions [20], and recording neuronal activities at >1.1 mm depth in the visual cortex of freely moving rats with a head-mounted miniature 3P microscope [21]. However, while the depth limit of 2PM has been extensively studied, systematic studies of depth limits of 3PM lack. Most of the current work that examines the depth limit of 3PM has employed theoretical simulations without experimental verifications [12,22]. The depth limit for 3PM has not yet been established theoretically or experimentally. New theoretical and experimental investigations are needed to fully understand the capability of 3PM and define the boundaries of high spatial resolution deep tissue imaging.
We follow the analytical model of Theer and Denk, where the spatiotemporal distribution of scattered and ballistic photons in a turbid media are analytically modeled to find a depth limit [8]. We extend this model to three-and four-photon (4P) excitation. Furthermore, we perform controlled experiments to quantitatively compare the depth limit of 2PM and 3PM and investigate the depth limit of 3PM in a variety of experimental conditions. Our theoretical and experimental data contradict published results from numerical simulation that showed loss of spatial resolution is the limiting factor in deep imaging with 3PM [22]. We demonstrate high spatial resolution imaging at depths of approximately 10 scattering mean free paths, which is nearly twice the transport mean free path (TMFP), the commonly acknowledged depth limit of imaging techniques using ballistic photons [23,24]. Our results unequivocally show that the depth limit of 3PM is significantly beyond what has been achieved in biological tissues so far. Our work will not only deepen the understanding of deep tissue imaging but will also motivate future technology development to reach the full potential of 3PM.

Theoretical analysis
To quantitatively describe the spatial distribution of fluorescence, we must find the distribution of excitation light inside a scattering medium. Ballistic photons are exponentially attenuated by tissue scattering and absorption, which are characterized by the scattering mean free path (l s ) and absorption length (l a ). The combined effect of absorption and scattering is represented by the effective attenuation length (EAL, l e ): The EAL is often used as the unit of depth in characterizing the optical imaging depth limit [8,15]. In our experiments, which mimic biological imaging conditions, scattering is the dominant cause of ballistic photon attenuation (i.e., l a ≫ l s and EAL ∼ l s ). This assumption holds for most biological tissues in vivo between 700 nm and 1700 nm (typical wavelength window for MPM) except for within the wavelength range of strong water absorption between approximately 1400 nm and 1550 nm.
The power of ballistic photons as a function of depth can be characterized by: where P 0 is the optical power at the surface of the media and z 0 is the distance traveled within the media. While absorption simply removes the excitation photons from the focal volume, scattering re-distributes the light intensity and increases the excitation outside of the focal volume. Multiphoton excited fluorescence signal from the focal volume (defined as the signal, S) is mostly generated by the ballistic photons. Therefore, the excitation power at the sample surface must exponentially increase as a function of depth to maintain sufficient power at the focus for detectable signal generation (Eq. (2)). As the imaging depth increases, the amount of fluorescence excitation outside the focal volume (defined as the out-of-focus background, B), which contains no high spatial resolution information, eventually overwhelms S. The signal-to-background ratio (SBR) characterizes the image contrast. The imaging depth limit can be defined by choosing a depth at which a required SBR value is reached. We chose the depth limit as the depth at which SBR = 1, which has been commonly used in the past [8,12,25].
To simplify the analytical modeling, the ballistic photons inside a sample can be formulated as an exponentially decaying Gaussian beam with the beam waist at the focus plane. An additional temporal component is added in the case of pulsed lasers to represent the temporal pulse shape. Assuming a Gaussian temporal distribution, the intensity of ballistic photons can be described as: where P is the normalized power through a transverse plane, z is the position along optical axis measured from the focus plane, ρ is the radial distance from the optical axis, ω(z) is the 1 e 2 radius of the beam, τ 0 is the 1 e 2 pulse half width, and α is the effective attenuation coefficient (1/l e ). When the beam is focused inside a scattering medium, it travels from a transparent medium to the absorbing-scattering environment of the sample. The light that travels inside the sample has two components: the ballistic photons that follow the Gaussian distribution (Eq. (3)) and the scattered photons that deviate from this path. The spatial distribution of light depends on many factors including the angular distribution of the scattered light. Scattering angular distribution is characterized with the anisotropy factor, g. In biological samples, scattering is mostly in the forward direction and the value of g is typically between of 0.8 and 1 [24].
The ballistic photons traveling inside the sample are governed by the focusing geometry and the properties of the sample. As the beam travels deeper inside the sample, the rays at larger convergence angles experience larger attenuation and cause the effective numerical aperture (NA) of the system to reduce with depth. The effective NA (NA eff ) in the sample at focus depth z can be approximated by using propagation laws of Gaussian beams as follows [8]: where NA = n sin θ, n is the refractive index of the medium, and θ is the far-field angular spread of the Gaussian beam. Using this approximation, we find the intensity distribution of the ballistic photons inside a scattering medium. We follow the model developed by Theer and Denk to find the intensity of the scattered light (I scattered ) [8]. A beam spread function is defined to analytically compute the mean and variance of the spatial and temporal distribution of the scattered light. The scattered intensity is then calculated by finding the convolution of the ballistic intensity with the beam spread function. Detailed description of the model is shown in Supplement 1.
Fluorescence produced by an m-photon process can be described as follows: The axial distribution of fluorescence (F m (z)) is found by integrating the radial and temporal components ( Fig. 1(a),(b))). To quantify the SBR given an axial intensity distribution, the signal is found by integrating the axial distribution of the fluorescence over increasing intervals around the focus depth, until a plateau in the value of S is reached (Fig. S1). We normalized the integration range to the Rayleigh length of the Gaussian beam as defined by As shown in Fig. S1, the signal reaches an asymptotical value for all focus depths. We chose a full integration range of 3z r and 5z r for 3P and 2P intensity distributions, respectively. The background value is determined by subtracting the signal from the total integrated fluorescence of the entire depth. The change in 2P and 3Ps S, B, and SBR with increasing focus depth is shown in Fig. S2.
An important factor in determining the background fluorescence is the fraction of the volume that is occupied by fluorophores since fluorescence is directly proportional to concentration [26]. This effect can be quantified by defining an inhomogeneity factor: Assuming the features of interest are larger than the focal volume, the background fluorescence can be scaled inversely with the inhomogeneity factor and the SBR can be calculated for samples with different inhomogeneity factors ( Fig. 1(c),(d))). Since fluorophore distribution in biological samples can vary, this representation of the depth limit allows us to adjust the values according to the sample properties. The depth limit will increase as the contribution of background decreases with increasing inhomogeneity, as shown in Fig. 1(d). Using the theoretical models described above, we have also investigated the effect of NA and pulse width of the excitation light on SBR ( Fig. S3 and S4).

Experimental setup
To verify the theoretical analysis above, we performed experimental measurements in tissue phantoms, which allow us to precisely control all the sample and imaging parameters. We embedded polystyrene beads in an agarose solution to create controlled absorbing and scattering tissue phantoms. We used 1 gram of high melting point agarose per 100 mL of water as our base medium so that we could stabilize the beads for resolution measurement while emulating biological samples with its absorption properties. To control the scattering and inhomogeneity factor we added a combination of 1-µm diameter fluorescent beads (Fluoresbrite Yellow Green Microspheres, 17154 Polysciences, Inc.) and 1-µm diameter non-fluorescent beads (Polybead Microspheres, 07310 Polysciences, Inc.). While keeping the agarose solution above its gelling point, the two types of beads were mixed with a vortex mixer, then mixed thoroughly with the agarose solution before pouring into the sample container to solidify.
Images were obtained using a commercial multiphoton microscope (Bergamo II system B242, Thorlabs Inc.). A high numerical aperture objective lens (XLPLN25XWMP2, Olympus, NA 1.05) was used for the experiments. The back aperture was approximately filled. The excitation source was a non-collinear optical parametric amplifier (NOPA, Spectra Physics) pumped by an amplifier (Spirit 1030-70, Spectra Physics). For 3PM and 2PM, the excitation wavelengths were 1280 nm and 920 nm, respectively. To reach the depth limit while maintaining low average power on the sample surface, low repetition rates of 333 kHz and 2 MHz were used for 3PM and 2PM, respectively. The full width at half maximum (FWHM) pulse width was measured to be 120 fs at 920 nm and 50 fs at 1280 nm, assuming sech 2 pulse intensity profile.
Images were collected at approximately 1 frame per second over a field-of-view (FOV) of 67 µm by 67 µm with 512 by 512 pixels. The point spread function (PSF) of both 2PM and 3PM images was characterized in shallow regions (<1 EAL) to verify diffraction-limited performance of the microscope. To find the EAL of the samples, z-stacks were taken. To represent the signal at the corresponding depth, the top 0.01% or 0.1% (depending on the inhomogeneity factor) of the brightest pixels of each image in the z-stack were averaged and normalized to the cubic or square (for 3P and 2P, respectively) of the power used. The signal value obtained by averaging the top 0.01% or 0.1% (depending on the inhomogeneity factor) of the brightest pixels of each image was used to calculate the SBR of each image (Fig. 2(c) and 4(c)). The background of each image was determined by averaging 90 percent of pixels with the lowest pixel counts. The slope of fluorescence signal versus depth in the semi-log plot was used to calculate the EAL [26]. All powers reported are values measured after the objective lens. For some samples, the EALs measured from the z-stacks were further verified with values measured with a spectrophotometer (Cary-5000-UV-Vis-NIR spectrometer, Agilent).

2PM and 3PM depth limit
We quantitatively characterized 2PM and 3PM images of tissue phantoms at various depths in order to experimentally determine the depth limits. A tissue phantom with a concentration of 9.4 × 10 6 beads per microliter was created using a combination of fluorescent and non-fluorescent beads. From the z-stack images, we measured the EALs to be 70 µm at 920 nm and 127 µm at 1280 nm ( Fig. 2(b)). The calculated scattering mean free path using Mie theory was 142 µm for wavelength of 1280 nm and 80 µm for wavelength of 920 nm. The specific bottle of nonfluorescent polystyrene beads used for this sample had a mean diameter value of 1.06 µm with a coefficient of variation of 3%. We attribute the difference between the measured and the calculated values to the variation in bead diameters. The Maximum powers of 74 mW at 1280 nm and 55 mW at 920 nm were used to obtain the deepest 3PM and 2PM images, respectively. Using the concentration of the bead solutions, an inhomogeneity factor of 4545 (∼0.02% fluorescent volume) was calculated for this sample. We verified the inhomogeneity factor experimentally by taking a fine z-stack of images spanning a portion of the sample depth and then counting the number of beads within the measurement volume. To quantitatively compare the experimentally measured and theoretically calculated SBR values at each depth, the theoretically calculated SBR at χ = 1 was scaled up by the inhomogeneity factor.
The SBRs of 2PM and 3PM images are similar in the shallow regions given the limited dynamic range of the microscope (Fig. 2(a)). At depths where 2PM images start to have noticeable background, 3PM has no measurable background generation. For this scattering sample, we reach an SBR of 1 at ∼560 µm (8 EALs) for 2PM, but with 3PM, we can penetrate deeper than 1270 µm, reaching beyond 10 EALs (Fig. 3(b)). As shown in Fig. 2(c), the SBRs measured in the sample in the deep regions are in good agreement with the theoretical predictions for both 2PM and 3PM. In the shallower regions, the theoretically calculated SBRs are too high to be experimentally verified with our microscope, particularly for 3PM. We note that the SBRs in these shallow regions are inconsequential for determining the depth limit of 2PM and 3PM. Individual beads are clearly resolvable even in the deepest images. To quantify the spatial resolution, we collected images of beads with a FOV of 6.7 µm by 6.7 µm with 256 by 256 pixels and at 0.3 µm z-steps in both the shallow and deep regions. The axial fluorescence intensity profiles are shown in Fig. 2(d). The FWHM are measured to be 1.6 µm for both 2PM and 3PM in the shallow region. The diffraction-limited resolution of the microscope for 2PM and 3PM is calculated to be 1.3 µm and 1.4 µm for NA = 1, respectively [27]. In the deeper regions the axial resolution degrades in high NA conditions due to attenuation of light at large convergence angles. The axial resolution in the deeper regions was measured to be 1.85 µm and 2.1 µm for 2PM and 3PM, respectively. The small degradation of the axial resolution (∼ 30% or less) in the deeper regions is expected due to the reduction of the effective NA when focusing in scattering samples. The lateral width of the beads was measured to have a FWHM of 0.75 µm and 0.8 µm in shallow regions and 0.82 µm and 0.79 µm in the deep regions for 3PM and 2PM, respectively. The diffraction-limited lateral resolution of our 3P imaging system is ∼ 0.4 µm (FWHM). The FWHM of the fluorescent beads in our images (∼ 0.8 µm) is noticeably smaller than the diameter of the beads (1.0 µm), which can be explained by the spherical geometry of the beads and possibly the non-uniform distribution of the fluorophores within the beads. The graininess of the PSF images was caused by photon shot noise. Figure 3 shows the decline of SBR in 3PM images of the sample. Images with SBR of above 1 were obtained at the depth of 10 EALs. We found close to diffraction-limited resolution near the depth limit (∼9 EAL deep in Fig. 2(d)), indicating that resolution degradation is not the limiting factor in reaching the depth limit. This finding contradicts a previous report that, based on numerical simulation, degradation of spatial resolution limits the imaging depth of 3PM [22].
With the anisotropy factor of the tissue phantoms calculated to be 0.82 using Mie scattering theory, the TMFP for the samples in our experimental conditions is [28]: Therefore, Fig. 3 shows that high spatial-resolution MPM can reach nearly twice the depth of the TMFP.

Inhomogeneity-dependent 3PM depth limit
To test the effects of inhomogeneity on the depth limit, we used samples with varying ratios of fluorescent to non-fluorescent beads. The tissue phantom with high inhomogeneity factor contained 9.9 × 10 6 beads per microliter and the tissue phantom with low inhomogeneity factor contained 9 × 10 6 beads per microliter. The scattering lengths were calculated using Mie theory to be 135 µm and 149 µm for the high inhomogeneity sample and the low inhomogeneity sample, respectively. We measured EALs to be 120 µm for the low inhomogeneity sample and 144 µm for the high inhomogeneity sample. The inhomogeneity of the samples differs by a factor of ∼13 (Fig. 4), which results in slightly less than 1 EAL difference in the 3PM depth limit. This result matches well with the theoretical calculation (Fig. 4(c)). Diffraction-limited resolution was obtained in both samples at ∼9 EALs deep in the sample.

Discussion and conclusion
MPM is a powerful tool for deep imaging due to reduction of the out-of-focus fluorescence and use of long wavelength photons that penetrate deeper into scattering tissue. While 2PM can reach deep within scattering tissue by using longer excitation wavelengths compared to one photon excitation, the optimal wavelengths for deep penetration in biological samples (i.e., ∼1300 nm and 1700nm) coincide with the 3P excitation wavelengths of the most commonly available fluorophores [5]. Additionally, the higher-order excitation in 3PM reduces the generated background fluorescence, producing high-SBR images beyond depths (normalized to the EALs) that are accessible to 2PM. While 4PM would further reduce background generation, we note that SBR is not the current experimental limitation for deep imaging with 3PM.
The TMFP is often stated as the fundamental limit for diffraction-limited optical microscopy since the ballistic photons are drastically attenuated by the time they reach this depth [23,24]. As we show with our experimental results, MPM can surpass this limit while maintaining a diffraction-limited resolution, even in densely labeled samples (e.g. with 3PM). Our results provide an intuitive understanding of three-dimensionally resolved optical microscopy as a competition between the in-focus signal and the out-of-focus background photons: the greater the reduction of the background photons, relative to the signal, the greater our chances of detecting the signal. While TMFP provides a length scale that describes the intensity distribution of the excitation beam, it does not characterize the nonlinearly generated signal and background. In MPM, the 'fundamental' depth limit is not caused by photons being able to reach a certain depth but rather the ability to distinguish the fluorescence photons that are generated at the focus from those that are from outside the focal volume, which is characterizable by the SBR.
The inhomogeneity factor is a key consideration in determining the depth limit since the amount of background generated by excitation light distribution is directly proportional to fluorophore distribution [26]. By linear fitting of the depth limit (normalized by EAL) as a function of logarithm of χ shown in Fig. 1(d) ranging from χ = 10 to χ = 10 5 , we found that the depth limit for 2PM and 3PM can be approximately expressed as: depth limit 3P EAL = 7.99 + 0.81 log 10 χ, depth limit 2P EAL = 3.67 + 1.29 log 10 χ.
Our theoretical results show that for approximately every 17 times increase in inhomogeneity, an additional EAL in depth will be accessible with 3PM. For 2PM the inhomogeneity needs to increase 6 times for each additional EAL of accessible depth. These results closely match our experimental findings. Equations (9 and 10) show that when the inhomogeneity factor is ∼ 10 9 , the depth limit for 2PM can eventually overtake that of 3PM. While such cross-over of 2PM and 3PM depth limits is theoretically interesting, such sparse labeling has little practical consequence. The inhomogeneity factor of 10 9 , which corresponds to having a single 1 µm fluorescent bead in a volume of 1 mm 3 in a tissue phantom, is most likely irrelevant to in vivo imaging. Therefore, within the practical range of in vivo imaging, the EAL-normalized depth limit of 3PM is larger than that of 2PM. By examining the differences in the depth limit ( Fig. 1(d)) it is clear that the advantage of 3PM over 2PM in penetration depth (normalized by the EAL) is more pronounced in densely labeled samples (i.e., samples with smaller inhomogeneity). For example, our simulations show that the normalized depth limit of 3PM for χ = 50 (similar to mouse cortex vasculature density) is approximately 9.4 EALs compared to depth limit of 5.9 EALs for 2PM, a difference of 3.5 EALs. For χ = 1000, the difference is reduced to approximately 2.9 EALs.
We plotted the axial fluorescence profile ( Fig. 1(b)), SBR as a function of depth ( Fig. 1(c)), and the imaging depth limit as a function of inhomogeneity ( Fig. 1(d)) using depth (z) normalized by the EAL. Such plots capture the effects of a variety of parameters on imaging depth limit. Intuitively, z/EAL represents the signal attenuation as a function of imaging depth (Eq. (3)) and is the most important parameter for defining the depth limit. On the other hand, the SBR also depends on the value of the EAL. Therefore, the coefficients in Eqs. (9 and 10) will depend on not only the imaging parameters (e.g., NA and pulse width) but also the EAL values. The EALs of biological samples can vary based on tissue properties. Fig. S5 shows that SBR increases with increasing values of EAL when plotted against the normalized depth, although the increase does not appear to be large enough to significantly change the depth limit (in units of z/EAL).
The deepest in vivo 3PM so far is at the depth of ∼6 EALs when imaging mouse brain vasculature [17]. Our results show that the current in vivo imaging depth achieved by 3PM is practically limited by the signal strength of 3P excitation. We used bright fluorescent beads in tissue phantoms to generate sufficient signal strength to reach the 3P depth limit. From the power levels at the surface of the sample, we estimated that the fluorescent beads we used are thousands of times brighter than typical in vivo imaging settings, e.g., imaging mouse brain vasculature. Therefore, brighter fluorophores together with schemes to enhance signal generation can lower excitation power requirements, making deeper regions more accessible [29,30]. Additionally, using lower repetition rate lasers [31,32] or adaptive excitation sources to reduce the effective repetition rate [33] can further improve the signal strength. Adaptive optics can also be used to enhance the 3P signal strength [34][35][36][37][38][39] by correcting the low order optical aberrations in biological samples. Combinations of some or all of these techniques can be employed to achieve the signal strength needed for reaching the 3P depth limit in vivo.
In summary, we extended the analytical model of Theer and Denk for 2PM to three-and four-photon excitation and performed controlled experiments to compare the depth limit of 2PM and 3PM. We also demonstrated the depth limit of 3PM in a variety of experimental conditions. We found that the analytical model is in good agreement with the experiments. Our experimental results show that loss of resolution (i.e., degradation of the point spread function) does not limit the imaging depth of 3PM. We show high spatial resolution at an imaging depth of > 10 scattering mean free paths, which is nearly twice the TMFP, the commonly acknowledged depth limit of imaging techniques using ballistic photons. We further studied the dependence of imaging depth limit on the staining inhomogeneity. Our results show that the depth limit of 3PM is significantly beyond what has been achieved in biological tissues so far. Our work will not only deepen the understanding of deep tissue imaging but will also motivate future technology development to reach the full potential of 3PM.