Two-photon excitation fluorescence microscopy (TPM) is a powerful tool for the in vivo three-dimensional (3D) imaging of cellular and subcellular structures and functions deep in turbid tissues [1]. Owing to its nonlinear excitation property, TPM provides compelling performance of near-diffraction-limited spatial resolution in relatively thick samples. However, the conventional TPM, i.e., Gaussian-focus TPM (Gauss-TPM), captures volumetric images by serially scanning a 3D space with a Gaussian focus, which significantly limits the imaging speed. Owing to its importance in capturing rapid biological events such as calcium transients in neurons, considerable effort has been devoted to improving the focus scanning rate [2,3], such as by applying resonant scanning [4], acoustic scanning [5–7], ultrasound lenses [8,9], electrotunable lenses [10,11], remote focusing [12,13], multifocal excitation [14–17], and multi-angle line scanning [18]. Other strategies, including temporal focusing [19,20] and multiplane imaging [21], have been developed to increase the volumetric imaging speed by capturing a two-dimensional (2D) image within a single exposure. Although these methods can partially alleviate the problem by improving the scanning speed, the implementation of layer-by-layer scanning to visualize 3D structures inherently limits the volumetric imaging speed. Moreover, the layer-by-layer scanning strategy repeatedly exposes the tissues above or beneath the focal plane to the excitation light, thereby aggravating photodamage and phototoxicity.

Sculpting the point spread function (PSF) of TPM to capture volumetric images may circumvent this problem [22]. Elongation of the PSF along the optical axis can be combined with raster scanning to form a projection image of a volume at the video rate [23]. This method obviates layer-by-layer scanning and advances applications such as functional imaging of neuron activity in sample volumes with sparse targets [24,25]. However, this method can only provide 2D projection images of the samples, whereas the axial location of fluorophores is not revealed. Meanwhile, the volumetric two-photon imaging of neurons using stereoscopy (vTwINS) [26] uses an elongated, V-shaped PSF to image a volume and encodes the depth location of a neuron into the separation distance between the paired projections in the resulting 2D image. This method records the neuron activity within a 45-μm-thick volume at a 30 Hz speed while preserving the depth location of the somas. However, the objective numerical aperture (NA) is partially sacrificed in vTwINS when creating a V-shaped PSF; consequently, the excitation efficiency and lateral resolution are degraded. Hence, this method is only suitable for imaging large somas, whereas the fine structures of neurons such as neurites cannot be inspected. Moreover, this method requires a sophisticated algorithm to demix fluorescent targets, thereby hindering its wide application in biological studies. Developing approaches that can increase the TPM volumetric imaging speed without sacrificing depth information or imaging resolution is vital to biomedical studies.

Herein, we present a simple yet versatile TPM volumetric imaging method based on paired-gradient excitation. In contrast to Gauss-TPM, the axial position of a fluorescent target is determined by the intensity ratio between the resulting paired images instead of by layer-by-layer scanning. This strategy increases the volumetric imaging rate by at least sixfold over Gauss-TPM and achieves a 0.63 μm axial localization precision, which is within the axial PSF of the Gauss-TPM; therefore, it is suitable for the rapid 3D imaging of sparsely labeled samples at subcellular resolution. In addition, our method significantly decreases the photodamage and phototoxicity caused by repeated exposure to out-of-focal plane excitation during axial scanning. Therefore, it is particularly suitable for in vivo imaging, such as the observation of macrophage phagocytosis dynamics.

The concept of the proposed scheme is straightforward. We used a phase-only spatial light modulator (SLM) conjugated to the back focal plane of the objective to generate a pair of axially elongated focal spots [Figs. 1(a) and 1(b)]. The intensity of the first focal spot decreased linearly along the optical axis with increasing depth, whereas the second focal spot exhibited an opposite intensity distribution. Subsequently, the specimen was scanned by the paired PSFs successively to obtain two z-projected images. Finally, the intensity ratio of the same component in the image pair was calculated, and a ratio-depth mapping function was used to determine the axial position of the structure [Fig. 1(a)]. We refer to this method as gradient TPM (Grad-TPM). In this system, an SLM is assembled into a Gauss-TPM system as an additional focus shaping module, which changes the Gaussian focus to a gradient focus with slight lateral resolution degradation (from 0.62±0.02μm to 0.76±0.08μm; see Fig. S1 in Ref. [27]). Switching between Grad-TPM and Gauss-TPM can be easily accomplished by loading specific patterns on the SLM. Throughout this study, a modest axial length of 12 μm was set for the gradient focal spot in Grad-TPM, resulting in a volumetric imaging speed improvement of at least sixfold, which depended on the z step size of Gauss-TPM. The z step size of our Gauss-TPM was set to 1 μm, considering the axial full width at half-maximum (FWHM) of the Gaussian PSF (1.2 μm) and the smallest repeatable incremental movement of the motor (0.8 μm).

To verify the feasibility of our approach, we compared volumetric images captured from Grad-TPM and Gauss-TPM in increasingly complex environments. First, we inspected the performance of our method by measuring 1-μm-diameter fluorescent beads embedded in agarose gel [Figs. 1(c)–1(e)]. A comparison of the two volumetric images shows that the lateral appearance of each bead in the Grad-TPM image was similar to the Gauss-TPM result. We further assessed the axial localization precision of our approach by comparing the axial center positions of beads from the Grad-TPM image pairs and the reference Gauss-TPM image stacks (1162 beads located at depths varying from 0 to 150 μm in the gel, obtained from 65 different imaging volumes). The results show that the axial localization precision, i.e., standard deviation of the localization errors, is 0.63 μm {Fig. 1(e) and Fig. S2 in Ref. [27]}, which is within the axial PSF of Gauss-TPM.

Next, we evaluated the performance of Grad-TPM on various biological samples. Volumetric imaging of neurons and neural networks is crucial for elucidating neural circuit functions. We imaged a fixed Thy1-GFP mouse brain tissue using Grad-TPM and Gauss-TPM, separately, and compared the 3D architecture of neural networks revealed by these two methods. A comparison between the Grad-TPM and Gauss-TPM images shows that the depth profiles derived along the axon were almost identical [Figs. 2(a) and 2(b)]. Moreover, a comparison of the corresponding sum intensity projection images shows similar intensity contrasts [Fig. 2(a) top], indicating that Grad-TPM is suitable for calcium imaging, which requires an accurate fluorescence intensity measurement of GCaMP indicators.

In addition to imaging fibrous structure in neural networks, cellular structure was also involved in the performance assessment of Grad-TPM. We measured the 3D distribution of microglia in CX3CR1-GFP brain slices [Figs. 2(c)–2(f)]. The average depth of each cell body calculated from Grad-TPM images was similar to the depth obtained from Gauss-TPM [the first row of Fig. 2(c) and Fig. 2(d)], and the depth variation along the microglia process was estimated well [the second row of Figs. 2(c) and 2(e)]. Besides, intensity contrasts of different subcellular structures were also well preserved in Grad-TPM [the third row of Fig. 2(c) and Fig. 2(f)]. Here, to capture multiple sparsely distributed microglia, we connected two axially adjacent volumes sequentially imaged by Grad-TPM [the first row of Fig. 2(c)]. The results show that the two volumetric images can be stitched seamlessly in the axial direction, demonstrating the potential of Grad-TPM for acquiring 3D images with greater axial extent.

After demonstrating that Grad-TPM can spatially localize objects in a 3D volume and accelerate volumetric imaging in phantom and fixed biological samples, we next investigated its photobleaching and phototoxicity in live cells. Photobleaching and phototoxicity induced by excitation light are two significant problems in the fluorescence microscopic imaging of biological samples from cells to organisms, particularly in living systems. To compare the photobleaching effects of Grad-TPM and Gauss-TPM, we conducted a time-lapse fluorescence imaging of live cells using the same volumetric imaging speed, the same initial fluorescence intensity, and the same frame rate. The detailed imaging conditions of the two methods are listed in Table. S1 in Ref. [27]. The results show that, after continuous imaging of HEK293 cells for 25 min, the fluorescence intensities of the cells imaged via Grad-TPM diminished slightly by 2.07±2.89% (number of cells, n=7), whereas in the Gauss-TPM, 70.85±8.93% (n=5) fluorescence was bleached [Fig. 3(a) and Visualization 1]. Even after 1.5 h of imaging, the reduction in fluorescence intensities measured from the Grad-TPM remained negligible [Fig. 3(a)].

Quantifying the phototoxic effects on live cells is less straightforward because the phototoxicity effect depends on many factors such as cell type, imaging conditions, intracellular distribution, and concentrations of fluorophores [28]. In this study, we assessed the phototoxicity effect by quantifying the rate of death of HepG2 cells during time-lapse imaging. To detect dead HepG2 cells, propidium iodide (PI), a classical dye that uniquely stains the nuclei of dead cells, was used. To locate live cells at the start of imaging, we simultaneously labeled the cells with DyLight 488 to highlight their membranes. For a relatively fair comparison, the excitation power was adjusted such that the peak fluorescence intensities from the same target were identical in the Grad-TPM and Gauss-TPM images, as the photobleaching experiment. Besides, the volumetric imaging speed and the frame rate used in Grad-TPM and Gauss-TPM were also set the same, respectively (see Table S1 in Ref. [27] for details). The results show that, in Grad-TPM imaging, none of the cells died even after 1.5 h of imaging, whereas the majority of cells died after 20 min of exposure to Gauss-TPM, suggesting that the phototoxicity was minimal during Grad-TPM imaging (Fig. S3 in Ref. [27] and Visualization 2). Although the reason for this is not fully understood, we believe that it is associated with the reduction in total exposure in Grad-TPM. Avoiding layer-by-layer scanning in Grad-TPM decreases the total irradiation dose on the samples. The photon dose of Grad-TPM for imaging the same volume was only one-third that of Gauss-TPM, despite the fact that the excitation power of Grad-TPM was approximately twice that of Gauss-TPM (Table S1 in Ref. [27]). The significantly lower photobleaching and phototoxicity effects of Grad-TPM on living biological samples render it more effective for obtaining reliable and reproducible quantitative data pertaining to biological processes compared with the standard 3D Gauss-TPM.

Owing to its rapid volumetric imaging speed and low phototoxicity, Grad-TPM is highly suitable for the longitudinal tracking of biological events in vivo. We used Grad-TPM to capture the 3D process of macrophages phagocytizing fluorescent beads at a speed of 6 volume/min [Figs. 3(b)–3(d) and Visualization 3]. We observed that individual beads agglomerated gradually to form a ring and were immobile subsequently, suggesting that the beads were phagocytized by the nearby macrophages and located inside the cell. To quantitatively analyze the bead motilities, we classified the beads into two categories, i.e., active beads and immobile beads, based on the criterion of whether the straight distance traversed by a bead from its starting point to the farthest point within 6 min exceeds 6 μm. The representative trajectories of these two types of beads are shown in Fig. 3(c). By calculating the mean velocity and plotting the mean displacement as a function of the square root of time, we discovered that the immobile beads exhibited constrained motility, whereas the active beads exhibited directed migration, although the velocities of these two beads did not differ significantly [29] [Fig. 3(d)]. These motility features indicate that the immobile beads might have been “locked” by the macrophages, whereas the active beads might have been transported into/in the macrophages. The results demonstrate the feasibility of using Grad-TPM to reveal the quantitative 3D dynamics and interactions of cells and macromolecules, which cannot be obtained via Bessel-beam-based volumetric imaging [23,24].

Limitations in the volumetric imaging speed and phototoxicity of Gauss-TPM considerably offset its advantages in the 3D imaging of living biological samples. Our methods solve this problem partially by modulating the excitation intensity along the z axis to generate a pair of gradient foci and integrating the depth information into the fluorescence intensity. This strategy can accelerate volumetric imaging by at least sixfold relative to Gauss-TPM while substantially mitigating photobleaching and phototoxicity by efficiently using excitation photons in the axial direction. Recently, Kumar et al. proposed a novel scanned oblique plane illumination (SOPi) microscopy method to provide rapid imaging with a speed up to 10 volumes per second [30]. It used a single objective to provide light-sheet scanning volumetric imaging. Compared with their relatively complicated system, our system can be easily implemented by adding an SLM conjugated to the objective back focal plane to the standard Gauss-TPM system without significantly modifying the excitation and emission paths.

In the experiments described herein, the axial length of the gradient focal spot in Grad-TPM was set to 12 μm, which is suitable for visualizing cultured live cells and brain tissue slices. The focal length can be easily adjusted from ∼4 to 20 μm by applying different phase patterns on the SLM (Fig. S4 in Ref. [27]). The maximum length was restricted by the SLM resolution, which determines the finest phase pattern that the SLM can generate, and the laser power, which determines the photon density of the extended excitation focus. Moreover, we discovered that excessively elongating the PSF magnified the axial intensity fluctuation of the gradient focus significantly, thereby decreasing the accuracy of axial localization. A modest length of 12 μm was optimal for our system setup. To perform deeper imaging, stitching in the axial direction, as shown in Fig. 2(c) (first row), is preferable over the further elongation of the PSF.

We note, however, several limitations or caveats in the current approach. First, we mapped the intensity ratio of the structures into their depth positions; therefore, the localization accuracy depended on the accuracy of the obtained intensity information. Optical scattering caused by refractive index differences in biological tissue will distort the linear gradient profile and accordingly lead to axial localization error. Optical aberrations caused by misalignment in the system and refractive index mismatch between the sample and the immersion media of the objective can further distort the gradient focus. Several approaches have been adopted to address these problems. Controlling the imaging depth in the superficial layer of biological tissue would effectively decrease the influence of optical scattering and optical aberration. Adaptive optics [31,32], which is used to compensate or reduce optical distortion, can restore the designed PSF in deep tissues. In addition, designing a more sophisticated axial intensity distribution instead of the simple linear distribution used in this study would enable the target components to be localized accurately. Second, the current scheme is only suitable for the visualization of sparsely labeled fluorescent samples, where only one fluorescent target is present in the axial direction. Otherwise, the axial locations of different fluorescence targets are averaged, and the targets cannot be distinguished well (Note S1 in Ref. [27]). This is a typical challenge in PSF axially extending methods for capturing volumetric images [23–26]. However, many biological structures, such as neurons and microglia, can be sparsely labeled in the imaging volume, as shown in this study (Fig. S5 in Ref. [27]). With our experience in testing kinds of biological samples, we think that samples with <15% overlapping ratio (the ratio of the area with axially overlapped structures within the 12 μm excitation range to the total fluorescent area, as shown in Fig. S5) are acceptable for the Grad-TPM. Please note that the condition of sparsity will decrease by increasing the focal length of the gradient foci. In addition, we are developing a stereoscopic Grad-TPM system that uses two pairs of V-shaped gradient foci for excitation. Therefore, the depth information is not only integrated into the fluorescence intensity but also encoded within the lateral shift patterns. In principle, this strategy can identify dense targets in a single volume. Third, as raster scanning was performed using galvanometer scanners, the volumetric imaging speed of our Grad-TPM system was only 30 volume/min. This rate is relatively modest, particularly when compared with the method based on an ultrasound lens [9]. On the other hand, since successive excitations of the pair of gradient foci were performed frame by frame rather than point by point (limited by the speed of the SLM), the axial localization accuracy of our method is easily affected by motions (such as cell migration in blood vessels) or fluorescence intensity fluctuations (such as calcium activity, noise, and photobleaching) in samples. Given 1 Hz frame rate, 0.76 μm lateral resolution, and 0.63 μm axial resolution of our system, axial localization information can be extracted without obvious artifacts only if lateral movement <0.76μm/s, axial movement <0.63μm/s, and intensity fluctuation <29.7%/s (see Note S2 in Ref. [27]). However, by incorporating Grad-TPM with fast scanning strategies, such as resonant scanning or polygon mirror scanning, and time-multiplexing excitation [16,17], the volumetric imaging speed may be improved to above 1800 volume/min, and the robustness to motions and intensity fluctuations will be significantly enhanced.

In summary, we proposed a new approach for capturing 3D images by encoding depth information into the fluorescence intensity. This method can substantially increase the volumetric imaging speed and decrease the photobleaching and phototoxicity of TPM with minimal resolution degradation. This concept, which uses the intensity ratio to encode the location information, can be incorporated into other nonlinear microscopies, including coherent anti-Stokes Raman scattering and stimulated Raman scattering, which are based on layer-by-layer scanning, to perform the 3D imaging of biological tissues.

A simplified diagram of the proposed Grad-TPM system is presented in Fig. 1(b). Pulsed light (920 nm, 80 MHz repetition rate) from a Ti:sapphire laser (Chameleon Ultra, Coherent), controlled by an electro-optic modulator (Model 350-80-LA-02, Conoptics), was expanded to a 1/e2 diameter of 30 mm before being reflected from the SLM (15° incident angle, 1920×1080pixels, PLUTO-NIR, Holoeye Photonics AG). A half-wave plate (GCL-060812, Daheng Optics) was placed before the SLM to orient the polarization angle of the incident laser to match accurately with that of the SLM. After being reflected by the SLM, the excitation light was raster scanned two-dimensionally by a pair of galvanometers (5 mm aperture; TS8203, Sunny Technology). Subsequently, an objective (1.15 NA, UAPON 40XW340, Olympus) was used to focus the scanning light on the specimen for eventual imaging. To ensure that the phase pattern from the SLM is stationary at the objective rear pupil plane, as the galvanometers were scanned, the SLM, paired galvanometers, and objective rear pupil plane were designed as mutually conjugate using three pairs of relay lenses (focal lengths fL3=150mm and fL4=75mm; fL5=fL6=85mm; fL7=80mm; fL8=200mm) (Fig. S6 in Ref. [27]). Furthermore, the relay lens pairs served as beam expanders or reducers to match the aperture of each element. In addition, a field stop located at the intermediate image plane between L3 and L4 was utilized to obstruct undesirable light, which was diffracted from the SLM [32]. In the detection path, the emission fluorescence signals were acquired by the same objective and then separated from the excitation light via a long-pass dichroic mirror that was placed immediately above the objective (T715LP, Chroma). Subsequently, a photon-counting photomultiplier (H7421-40, Hamamatsu) was used to detect the fluorescent signals that were refocused by a collection lens (f=75mm) and spectrally filtered by appropriate bandpass filters (MF525-39 or MF620-52, Thorlabs). Scanning and data acquisition were controlled with a commercial data acquisition system (PCIe-6110, National Instrument) in combination with a custom-designed LabVIEW (National Instruments) software. For the extended Grad-TPM 3D imaging that connects multiple axially adjacent volumes or Gauss-TPM volumetric imaging based on layer-by-layer scanning, the objective was axially translated using an actuator (KMTS25E, Thorlabs).

In this study, the gradient excitation focus pair was designed to possess axial intensity distributions that were not only linearly varying but were also complementary. To engineer the designed gradient focal spot, the phase of the incident wavefront at the pupil plane (i.e., pupil phase) was manipulated by an SLM. Conventional pupil phase calculation methods based on the paraxial theory of scalar diffraction, such as Gerchberg–Saxton phase retrieval [33] and the extended Nijboer–Zernike theory [34], perform well in low-NA optical systems but imprecisely in high-NA configurations. Thus, we developed an algorithm that calculates the pupil phase of gradient focuses for high-NA systems. The kernel of the algorithm includes three components: ring segmentation to facilitate the establishment of the pupil phase function; the Richards–Wolf vector diffraction theory, which correlates the intensity distribution of the focal spot and the pupil phase function; and the genetic algorithm (GA) for searching the optimal pupil phase pattern that yields the desirable gradient focus (see Fig. S7 in Ref. [27]).

An ideal two-photon microscope can be regarded as a rotationally symmetrical system. Therefore, the pupil phase function expressed in polar coordinates can be simplified considerably by using only radial coordinates in the calculation. We restate the pupil as a series of area-equivalent concentric rings, each having a phase φk(k=1,2,…) to form the pupil phase function, i.e., φk(k=1,2,…,40) [35]. The intuition underlying ring segmentation is to manipulate the relative phases of light beams in each concentric ring and recombine their light fields at focal zone to form a predesigned intensity distribution. By increasing the number of rings, more precise pupil phase functions can be obtained. Forty rings were adopted in our system. They cover a circular area with a diameter of 1000 pixels together on the SLM, approximating to the short side of the SLM (1080 pixels). We found 40 rings can provide enough accuracy to generate a 12-μm-long linear gradient focus pair, as well as save computational resources.

The Richards–Wolf vector diffraction theory [36,37], a generalization of the Debye integral representation, can precisely describe the electromagnetic field in the image space of a high-NA optical system, for which the scalar diffraction theory is deficient. Based on this theory, we established a relationship between the pupil phase function and the focal field and, hence, the focal intensity distribution. To simplify the computation, only the axial intensity in the range 25 μm above and below the objective focal plane was considered in the calculation. We assume that the illumination is uniform and linearly polarized; therefore, the axial excitation intensity distribution of the focal spot can be calculated as I(z)=|E(z)|2 using the Richards–Wolf theory, where E(z) is the axial field distribution expressed as [38] E(z)=∫0αP(θ)(cosθ)12sinθ(1+cosθ)e−ikzcosθdθ.(1)P(θ)=φk(θk−1<θ<θk) is the pupil phase function based on ring segmentation, where θ is the angle of convergence of the objective lens, namely the aperture angle, and θk is the upper-bound aperture angle of the kth ring. α is the maximum θ, which is determined by the NA of the objective lens and the medium refractive index, n, based on α=arcsin(NA/n); k=2π/λ is the wavenumber, where λ is the wavelength. In our scheme, NA=1.0, n=1.33, and λ=920nm. The calculation of I(z) from P(θ) is provided in supplement Code 1, Ref. [39].

GAs are stochastic global search and optimization methods that mimic the metaphor of natural biological evolution [40]. GAs are suitable for large-scale optimization problems and are, therefore, attractive for optimizing phases in the focusing task [41]. In this study, we adopted a GA provided by Chipperfield et al. [42]. to determine the pupil phases of the 40 subregions partitioned via ring segmentation. As shown in Fig. S7 in Ref. [27], our GA first generated an initial population of 8000 individuals. Each individual is a pupil phase function with 40 chromosomes, which correspond to the phases of the 40 rings. Subsequently, the I(z) of each individual was calculated based on the Richards–Wolf theory, and their similarity to the target profile It(z) was measured using an objective function, Fobj=∑z[I(z)−It(z)]2, to assess the fitness of these individuals for the subsequent selection. Both I(z) and It(z) were normalized to their maxima. Next, the GA generated a new population for the next round of iteration through selection, crossover, and mutation, all of which were performed via default functions provided by the toolbox. The size of the population decreased to 400 when the number of generations reached 10 and was further reduced to 50 when it reached 50. The large initial population size prevented convergence to a local maximum, and the latter reduced size accelerated the algorithm. Eventually, the iteration stopped when the number of generations reached 4000, and the individual with the minimum Fobj in the last population was regarded as the global best satisfactory individual.

The gradient focus generation algorithm can converge to a globally optimized pupil phase function, although the process is time-consuming. Nevertheless, the imaging speed will not be restricted because the processes of pupil phase optimization and Grad-TPM imaging are completely independent. Finally, using the method above, a pair of gradient foci encompassing a 12 μm depth was obtained in the simulation and confirmed experimentally (Figs. S8 and S9, Table S2 in Ref. [27]).

To perform Grad-TPM imaging, the specimen was sequentially scanned using the gradient focus pair, Focus 1 and Focus 2, which generated an image pair, Im1 and Im2, containing depth information in their intensity. More specifically, Focus 1 exhibited an intensity that decreased linearly with increasing depth, whereas Focus 2 exhibited an opposite intensity distribution. Therefore, the axial location information was linearly encoded into the intensity of the gradient focus. Because the intensity of two-photon excitation fluorescence is proportional to the square of the excitation intensity, the axial location z of a fluorophore is further encoded into its emission intensity and, hence, can be decoded from the obtained image pair, as follows: z=f(R),R=Im2Im1+Im2,(2)where R is the ratio of excitation intensity of Im2 to the total. Function f(R) is a ratio-depth mapping function, which maps R to z. In principle, f(R) is a linear function as shown in Fig. 1(a). But actually it is not exactly linear, since the axial intensity of the generated gradient excitation focus is not ideally linear (Fig. S9b in Ref. [27]). To achieve a more precise localization, we corrected f(R) using the experimentally measured intensity profiles of the focuses. By applying function f(R), a z matrix corresponding to the imaged volume can be obtained.

Because the axial location information is deciphered from the intensity ratio, the noise contained in Im1 and Im2 significantly decreases the axial localization precision. To circumvent this problem, a three-step strategy was adopted, as shown in Fig. S10 in Ref. [27]. First, the background was eliminated, and only the pixels that contained objective fluorophores were used for analysis. This was performed by first normalizing Im1 and Im2 to their maximum intensities and then merging them to form a homogeneous excitation image Im that reflects the true intensity of fluorophores in the specimen, as follows: Im=Im1+Im2+2Im1·Im2.(3)Subsequently, a binary mask of Im, which was created using a threshold determined using Otsu’s method, was applied to Im to isolate the objective fluorophores from the background. Next, fluorophores outside the designed axial range of measurement were identified and excluded. As shown in Fig. S9b in Ref. [27], the gradient focus actually exhibits intensity outside the designed axial range, which may also excite some fluorophores, and their axial locations will be incorrectly estimated. By defining the upper and lower limits of the intensity ratio R, the fluorophores located above or below the designed range were ascribed to the upper and lower boundaries of the designed axial range, respectively. The detailed process is shown in Fig. S10 in Ref. [27]. Subsequently, we applied a moving average filter with a window size of 3×3pixels to smooth the z matrix and decrease the random noise.

Eventually, to visualize the volumetric imaging result of Grad-TPM, we created a color-coded image that preserves both the spatial location and the intensity contrast information of the specimen by coding the depth, z, with colors and the intensity, Im, with color saturation (Fig. S10f in Ref. [27]).

Gauss-TPM was used to provide ground-truth z values for comparison with the Grad-TPM in Figs. 1 and 2. For fluorescent beads in Fig. 1, we fitted the axial emission intensity profile of each bead with Gaussian fitting to determine their axial center (where the maximum fluorescence intensity was located). For biological samples in Fig. 2, we used a weighted summation method to locate the axial centers of the imaged features (zc): zc=∑z=112zρ(z)∑z=112ρ(z),(4)where z is the depth of each frame and ρ(z) is the fluorescence intensity obtained at the depth of z.

To measure the intensity distribution of the gradient focus and evaluate the lateral resolution and axial localization precision of Grad-TPM, fluorescent beads embedded in agarose gel were used. To prepare the sample, 1–10 μL suspension of 1-μm- or 0.1-μm-diameter yellow–green fluorescent beads (1:100 dilution, F8823 or F8803, Invitrogen) was mixed with 1 mL of 1% agarose solution. The mixture was vortexed vigorously and then deposited on a No. 1.5 glass bottom dish (P35G-1.5-14-C, Matek). The gel was allowed to solidify for several minutes and then imaged immediately.

To evaluate the performance of Grad-TPM on biological samples, fixed brain slices from Thy1-GFP and CX3CR1-GFP transgenic mice (∼4 weeks old) were prepared for imaging. The mice were maintained in a specific-pathogen-free animal facility of the Shenzhen Institutes of Advanced Technology (SIAT), Chinese Academy of Sciences (CAS). Thy1-GFP mice express EGFP primarily in mossy fibers in the internal granule layer of the cerebellum. CX3CR1-GFP mice express EGFP in the brain microglia, monocytes, dendritic cells, and NK cells.

To prepare the brain slices, the mice were deeply anesthetized with a mixture of 2% α-chloralose (30037517, Sinopharm) and 10% urethane (8 mL/kg, 30191228, Sinopharm), via intraperitoneal injection. Transcranial perfusion with PBS and 4% (mass-to-volume ratio) paraformaldehyde (PFA) in PBS was performed. The mice were sacrificed using this operation. Subsequently, the mouse brain was excised and fixed with 4% PFA at 4°C overnight. Finally, 500-μm-thick coronal slices were sectioned freehand using a brain matrix (68713, RWD Life Science). The brain slice was then placed on a glass slide and covered with a microscope coverslip for immediate imaging.

All experiments were performed in compliance with protocols that had been approved by the Guangdong Provincial Animal Care and Use Committee as well as guidelines of the Animal Experimentation Ethics Committee of SIAT, CAS.

Three cell lines, HEK293, HepG2, and RAW264, were used in this study. All cells were cultured in high-glucose DMEM (SH30022.01, HyClone) supplemented with 10% fetal bovine serum (FBS, 2023-02, Gibco), 50 U/mL penicillin, and 50 μg/mL streptomycin (SH40003.01, HyClone) in an incubator at 37°C with 5% CO2.

For photobleaching evaluation, HEK293 cells were transfected with the pCAG-EGFP construct using Lipofectamine 2000 (#11668-027, Invitrogen). pCAG-EGFP DNA (200 ng) was used in every 34 mm glass bottom dish (D35C4-20-1.5-N, In Vitro Scientific). The imaging buffer was phenol red-free DMEM/F12 (Gibco, 11039-021) containing 10% FBS.

For phototoxicity assessment, HepG2 cells were labeled with PI (P21493, Thermo Fisher) and DyLight 488 (DL-1174, Vector Lab). For staining, cells at a density of 70%–80% were incubated in a mixture of PI (1:3000 dilution, 500 nmol/L) and DyLight-488-labeled Lycopersicon esculentum (Tomato) lectin (10 μg/mL) for 10 min.

To observe the phagocytosis of macrophages, RAW264 cells mixed with fluorescent beads were imaged. Cells at a density of 70%–80% were passaged into two new dishes and incubated continuously for 4–5 h. Subsequently, 10 μL suspension of 1-μm-diameter yellow–green fluorescent beads (F8823, Invitrogen) was added to the incubation solution and blended with the cells. Finally, after incubating the mixture of cells and beads for 5 min, the sample was imaged immediately.

We used the Richards–Wolf vector diffraction theory to calculate the electric field intensity in the x−z plane. We measured the FWHM of the gradient focus at different z positions, and the simulation results show that the lateral resolution of the Grad-TPM is approximately 0.72 μm theoretically. Furthermore, we acquired the experimental lateral resolution by imaging 100-nm-diameter yellow–green fluorescence beads in the gel. The statistical result shows that the experimental lateral resolution of Grad-TPM is 0.76±0.08μm (Fig. S1 in Ref. [27]).

To assess the performance of Grad-TPM, particularly in terms of axial localization precision and preservation of contrast information, Gauss-TPM was used as a standard 3D imaging method for comparison. For all imaging experiments, a 920 nm laser was used for excitation, and a 505–545 nm bandpass filter was used for the spectral filtering of signals before detection, except for imaging DyLight 488 labeled membrane of HepG2 cells during phototoxicity assessment. DyLight 488 was excited by a 960 nm laser, and the emission fluorescence signals were spectrally filtered using a 594–646 nm bandpass filter. For photobleaching and phototoxicity assessment, the values of the initial fluorescence intensity, pixel dwell time, field of view, and pixel size were maintained for Grad-TPM and Gauss-TPM imaging. To obtain a similar fluorescence intensity between Grad-TPM and Gauss-TPM images, the excitation laser power of Grad-TPM was set to 2 to 3 times that of Gauss-TPM. The data acquisition parameters for each experiment are provided comprehensively in Table S1 in Ref. [27].

Classical motility analysis methods were used to quantify the dynamics of fluorescent beads during macrophage phagocytosis. Beads that were trapped by cells were used for motility analysis. We tracked the beads by recording their 3D centroid coordinates in each time-lapse image to obtain their trajectories and then calculated the motility parameters, such as the velocity and displacement. All analyses and calculations were performed using custom MATLAB (MathWorks) programs.

For data with a normal distribution, mean and standard deviation were used for estimations, while for other data with unknown distribution, boxplots that show the interquartile ranges were provided. No statistical method was used to predetermine the sample size. Data collection and analysis were not blinded to the experimental conditions. None of the data were excluded from the analysis.

Acknowledgment. We thank Dr. Hari Shroff and Dr. Yicong Wu for their helpful discussion and suggestion.