Adaptive optical microscopy via virtual-imaging-assisted wavefront sensing for high-resolution tissue imaging

Introduction Adaptive optics (AO) is widely employed in optical microscopy to achieve ideal spatial resolution in biological tissues [1, 2]. Correction of sample-induced aberrations by AO restores not only the resolution by recovering the aberrated wavefront but also the signal intensity by focusing the light more efficiently. Measuring the aberrated wavefront is a necessary step for AO, which can be generally categorized into direct wavefront sensing or indirect wavefront sensing [3, 4]. Direct wavefront sensing methods use a dedicated wavefront sensor, such as the Shack-Hartmann wavefront sensor, to measure the distorted wavefront originating from a guide star inside biological tissues [5–7], and these methods are usually applied to transparent or weakly scattering samples. In contrast, indirect wavefront sensing methods determine the aberrated wavefront from a sequence of images or the variation in the signal intensity instead of a wavefront sensor and thus have simpler hardware implementation and are preferable for scattering samples. Abstract


Introduction
Adaptive optics (AO) is widely employed in optical microscopy to achieve ideal spatial resolution in biological tissues [1,2]. Correction of sample-induced aberrations by AO restores not only the resolution by recovering the aberrated wavefront but also the signal intensity by focusing the light more efficiently. Measuring the aberrated wavefront is a necessary step for AO, which can be generally categorized into direct wavefront sensing or indirect wavefront sensing [3,4]. Direct wavefront sensing methods use a dedicated wavefront sensor, such as the Shack-Hartmann wavefront sensor, to measure the distorted wavefront originating from a guide star inside biological tissues [5][6][7], and these methods are usually applied to transparent or weakly scattering samples. In contrast, indirect wavefront sensing methods determine the aberrated wavefront from a sequence of images or the variation in the signal intensity instead of a wavefront sensor and thus have simpler hardware implementation and are preferable for scattering samples.
Among indirect methods, modal wavefront sensing methods use a quality metric, such as brightness, of the recorded image as a clue to retrieve the aberrated wavefront. By intentionally distorting the wavefront with a series of orthogonal aberration modes, such as Zernike polynomials, the amplitudes of these aberration modes of the aberrated wavefront can be calculated according to the variation in the resulting image metric [8][9][10]. In addition to modal wavefront sensing, the aberrated wavefront can also be obtained in a zonal manner, in which the entire aberrated wavefront is divided into segments based on pupil segmentation and reconstructed from the gradient of each constituent wavefront segment [11]. A wavefront segment of the aberrated wavefront at the back pupil of a microscope objective corresponds to a focusing beamlet that is deflected from the ideal focus of the objective, and the gradient of this wavefront segment can be calculated from the lateral displacement between the beamlet focus and the ideal focus. By acquiring an image under full-pupil illumination as a reference image, the lateral displacement of the beamlet focus can be reflected as the shift in the image acquired with that beamlet (shifted image) relative to the reference image [11]. Therefore, the gradient of a wavefront segment can be calculated from the shift between the shifted image and the reference image. However, for a wide variety of samples, since the beamlet focus is substantially elongated due to the much reduced effective numerical aperture (NA), the shifted image captures sample structures beyond the excitation volume of full-pupil illumination. As these additional structures dominate the shifted image, deriving the underlying shift between the shifted image and the selected reference is very difficult. To enable aberration measurement via image shift measurement in these samples, one can introduce fluorescence beads into the biological samples, and then, the aberrated wavefront can be measured using the sparsely distributed exogenous beads rather than the intrinsic biological structures [12,13]. Instead of image shift measurement, the gradient of each wavefront segment can also be obtained by applying a set of gradients to this wavefront segment and retrieving the gradient that leads to the maximum fluorescence intensity [14]. This approach is not limited by the sample structure, and the gradient measurement for different wavefront segments can be parallelized by using a high-speed wavefront modulator, such as a deformable mirror (DM), to implement a frequency multiplexing scheme [15]. However, since the gradient of a wavefront segment is determined based on the maximum signal intensity during wavefront modulations, the signal intensity without the intended modulations needs to be constant.
Here, we report an indirect wavefront sensing method that leverages a virtual imaging scheme to measure aberrations using intrinsic structures by image shift. To derive the lateral displacement of each beamlet focus from the resulting shifted image, the virtual imaging scheme computationally constructs an ideal beamlet focus, which corresponds to the same wavefront segment as the actual beamlet and is free from aberrations, based on the focusing model of the objective. The sample structures are virtually imaged via computations using the constructed ideal focus according to the imaging model. Since the ideal beamlet has the same effective NA as the actual beamlet, most of the content in the shifted image corresponds to a shifted counterpart in the virtually imaged reference image, which is crucial for measurement of the underlying shift. The underlying shift is then measured with a structural-similarity-based shift measurement algorithm and used to derive the corrective wavefront. We implement this method in a two-photon fluorescence microscope with a commonly used liquid crystal spatial light modulator (SLM) as a wavefront shaper and demonstrate the effectiveness of our method in a variety of samples, from dense fluorescent beads to fixed biological tissues and living animals. Since SLMs are also widely used in resolution enhancement techniques [16][17][18][19][20][21], our method can be directly integrated with these techniques to obtain higher resolution in biological tissues without any hardware modifications, and we present AO-incorporated subtractive imaging as an example. The robustness of our method to signal variation is further demonstrated by simulations and aberration measurement on neurons exhibiting spontaneous activity in a living larval zebrafish.

Principle of wavefront measurement via virtual imaging
To measure the aberrated wavefront induced by biological tissues, the wavefront at the back pupil of the microscope objective is split into an array of segments (Methods). The gradient of each wavefront segment is measured first to obtain the entire aberrated wavefront. As shown in Fig. 1a, according to the focusing model of the objective, if a wavefront segment is tilted with gradient ⇀ g , then the focus of the corresponding beamlet will be laterally displaced relative to the ideal focus of the objective along vector − → d , which can be calculated as where f is the focal length of the objective and λ is the wavelength. In point-scanning-based microscopy, where the image is formed by scanning the focus across the sample, the lateral displacement of the focus results in the same displacement of the field of view (FOV), which is consequently presented as a shift of the resulting image. Therefore, by acquiring an image with the beamlet corresponding to a wavefront segment (shifted image), the gradient of this wavefront segment can be estimated from the shift ⇀ s between the shifted image and a reference image acquired with the corresponding ideal beamlet, which is free from aberrations, as ⇀ g = c ⇀ s , where c = 2πl/(fλ) and l is the pixel size of the image. A schematic of measuring the gradient of a representative wavefront segment in biological tissues is depicted in Fig. 1b. A plane with the feature of interest, such as a neuron, is chosen as the central imaging plane for aberration measurement (e.g., the z = 0 plane in Fig. 1b-I). By illuminating a small zone at the back pupil of the objective corresponding to the wavefront segment and laterally scanning the resulting beamlet, which is deflected from the ideal focus by the sampled-induced aberration, across the sample, the structures excited at the beamlet focus are recorded in the shifted image ( Fig. 1b-II). Ideally, the gradient of the wavefront segment could be derived from the shift between the shifted image and a reference image acquired with the corresponding ideal beamlet. As the ideal beamlet is not available in practice, the conventional approach [11] acquires an image under full-pupil illumination as the reference image. Based on this reference image, the underlying shift can be measured from the shifted image when the sample structures are sparsely distributed (Fig. S1a). However, since the beamlet focus is greatly elongated due to the low effective NA, the shifted image can record structures beyond the excitation volume of full-pupil illumination. For a variety of samples, these additional structures are dominant in the shifted image, so the underlying shift cannot be derived from the selected reference image acquired under full-pupil illumination (Fig. S1b). Therefore, a reference image containing structures similar to the shifted image is essentially needed for accurate shift measurement. To this end, we leverage a computational scheme, termed virtual imaging (Methods), to construct the reference image. Prior to virtual imaging, the sample structures surrounding the central imaging plane are imaged under full-pupil illumination and recorded as a three-dimensional stack of images ( Fig. 1b-III). Next, the ideal beamlet focus, which corresponds to the same wavefront segment but without aberration, is calculated based on the focusing model of the objective and used to virtually image the recorded structures through computations. The computed image is then taken as the reference image for shift measurement (Fig. 1b-IV). Because the effective NA of the ideal beamlet is the same as that of the actual beamlet, most of the contents in the shifted image have a shifted counterpart in the constructed reference image, which provides the basis for accurate shift measurement. After the reference image construction, we use a structural-similarity-based shift measurement algorithm (Methods) to obtain the shift between the shifted image and the reference image ( Fig. 1b-V). This algorithm places the shifted image at every position within the reference image and measures the structural similarity between the shifted image and the overlapped patch in the reference image each time. The measured results can be viewed as a structural similarity map, with which the underlying shift can be estimated from the position with the largest similarity. The gradient of the wavefront segment can thus be obtained from the calculated shift.
With the gradient of each wavefront segment of the aberrated wavefront measured, we can obtain the corrective wavefront by applying wavefront reconstruction algorithms (Methods) used for the Shack-Hartmann wavefront sensor [22,23]. As shown in Fig. 1c, two corrective wavefronts are reconstructed by zonal and modal reconstruction algorithms. To obtain a better correction performance, each of these corrective wavefronts is applied for aberration correction. The corrective wavefront that gives a higher signal intensity under the same input power for the objective is then determined as the final corrective wavefront.

Wavefront measurement with dense fluorescent beads
We implemented our method in a homebuilt two-photon microscope (Methods), employing an SLM to achieve zonal illumination and wavefront correction. To verify the ability of the proposed method to measure an aberrated wavefront with densely distributed structures, we performed AO correction for a dense aggregate of 2 μm fluorescent beads embedded in agarose with applied artificial aberration. To simulate the aberrations induced by biological tissues, we generated the artificial aberrations from the experimentally measured aberrations in biological tissues by adding perturbation to the Zernike coefficients.
We first captured a sequence of shifted images for different wavefront segments under zonal illumination (Fig. 2a). Compared with the image acquired under full-pupil illumination (indicated by the cyan square in Fig. 2a), the shifted images are dominated by additional structures beyond the excitation volume of full-pupil illumination. To derive the underlying shifts of these images, we acquired an image volume of the beads surrounding the original imaging plane (Fig. 2b) and constructed reference images based on the acquired volume through the virtual imaging scheme. For each shifted and reference image pair, the corresponding similarity map was measured, and the underlying shift was estimated from the position with the largest similarity (Fig. 2c). As shown in Fig. 2d, the image segment that corresponds to the estimated shift of the reference image matches the shifted image, indicating that the underlying shift was successfully measured. Before AO correction, the beads in the image were dim and elongated in the axial direction ( Fig. 2e (No AO)). When the AO correction was performed with the Virtual-imaging-assisted AO enables aberration correction for dense 2 μm fluorescent beads. a Shifted images acquired under zonal illumination and image acquired under full-pupil illumination. The entire wavefront is split into an array of segments (red dashed squares), and the acquired shifted images are arranged by the corresponding wavefront segment. The image marked by the cyan square was acquired under full-pupil illumination with the same FOV as the shifted images, and it was used as the reference image in the conventional method. b Maximum intensity projection of the captured image volume, which serves as the sample in virtual imaging. The yellow dashed box in the lateral view marks the imaging FOV during zonal illumination. The yellow dashed line in the axial view indicates the objective focal plane during the aberration measurement process. c Maps of structural similarity with respect to the shift. In each map, the shift with the largest similarity is marked by the red dot. d Shifted image, reference image and similarity map for a representative wavefront segment marked by the yellow asterisks in (a) and (c). The red box indicates the patch of the reference image exhibiting the largest similarity with the shifted image. e Maximum intensity projection of the image volume without AO correction (No AO), with AO correction using the conventional method (Conv) and our method (AO), and under the initial condition where no artificial aberration is applied (Init). The insets in the images without AO correction and with AO correction using the conventional method are digitally enhanced 3-fold and 8-fold, respectively, to increase the visibility. f Corrective wavefront. The dashed squares indicate the layout of the wavefront segments, where each wavefront segment corresponds to a unit square. An independent mask approach was used in the measurement process. g Signal profiles along the cyan, brown, magenta and orange lines in e. Scale bars, 5 μm conventional method, which uses the image acquired under full-pupil illumination as the reference, the correction led to a deteriorated imaging performance, as evidenced by a 1/3-fold decrease in the peak signal intensity and severe image distortion ( Fig. 2e (Conv)). In contrast, after applying our method (the corrective wavefront is shown in Fig. 2f ), both the signal intensity and shape of the beads were almost restored to the original state when no artificial aberration was applied ( Fig. 2e (AO, Init)). According to the intensity profiles shown in Fig. 2g, AO correction using our method reduced the full width at half maximum of the bead from 4.2 μm to 2.6 μm and increased the signal intensity by nearly 3-fold. The imaging performance after AO correction approaches that when no artificial aberration was applied, which verifies that our method can be applied to densely distributed structures.

Aberration correction in biological tissues
To test the effectiveness of our method in biological tissues, we imaged neuronal structures expressing the red fluorescent protein mRuby3 in a fixed mouse brain tissue section. We measured the aberration with the neuronal structures ( Fig. 3a, b) axially centered 250 μm below the tissue surface and observed substantial improvements in both image brightness and contrast after AO correction (Fig. 3c). As shown in Fig. 3d, AO correction increased the fluorescence intensity of a soma (the pink line in Fig. 3c) and a dendrite (the cyan line in Fig. 3c) by approximately 2.5-fold and 4-fold, respectively, and allowed more neuronal structures to be observed. As the axial resolution is more sensitive to the presence of aberrations, the resolution improvement brought by AO correction can be readily recognized from the axial views in Fig. 3c, e, where more dendrites can be clearly distinguished after AO correction. With the signal gain and the resolution improvement provided by AO correction, dendritic spines that are unresolvable without AO correction (the yellow arrows in Fig. 3f ) can also be identified (Fig. 3f ).
To explore the capability of our method for living scattering samples, we applied our method to imaging the tumor microenvironment in a living mouse with enhanced green fluorescent protein (EGFP) expressed throughout the entire body, excluding erythrocytes and hair. mCerulean-B16 tumor cells were injected between the fascia and dermis of the rear skin of the EGFP mouse, and the tumor region was imaged through a skin-fold window chamber implanted on the back of the EGFP mouse. Prior to imaging, elastin fibers of the arterial vessels were labeled by Alexa Fluor 633 via tail vein injection so that the tumor cells, host-derived EGFP-expressing cells and arterial vessels could be distinguished by three different fluorescence spectra (Fig. 4a).
We measured the aberration with the microvessels axially centered 92 μm below the tissue surface (Fig. 4b). As shown in Fig. 4c, enhanced fluorescence signals were observed after AO correction in all three channels, where the signal intensities of a host cell (pink line), a tumor cell (orange line) and a microvessel (cyan line) were enhanced approximately 1.7-fold, 2-fold and 3-fold, respectively (Fig. 4d). The signal enhancement would be helpful for better distinguishing cell types and quantifying the number of cells. With AO correction, the resulting signal gain and resolution improvement allowed more microvessels to be visualized (Fig. 4e) and individual microvessels to be resolved from a blurred cluster (Fig. 4f ). Thus, AO correction enables observation of the detailed network structure of microvessels, which would benefit the monitoring of the morphological development of tumor microvessel network.

AO-incorporated subtractive imaging
Resolution enhancement techniques have been proposed to pursue higher resolution in biological tissues, and many of them employ SLMs to implement the desired wavefront modulations [16][17][18][19][20][21]. Since only a single SLM is required in our AO correction method, our method can be readily incorporated with these techniques without extra hardware complexity. Here, as an example, we demonstrate the application of our method with subtractive imaging [21,24,25], which enhances the resolution and contrast of the conventional imaging modality by image subtraction. In brief, subtractive imaging employs a conventional focal spot and a donut-shaped focal spot, which can be shaped through wavefront modulation, to probe the sample (Methods). By leveraging the different feature sizes for the focal spots, an image with enhanced resolution and contrast can be constructed via weighted subtraction between images acquired with the different focal spots.
We imaged a Pinus pollen grain (Pine Mature Pollen, Carolina Biological Supply Company), as shown in Fig. 5a. Before applying subtractive imaging, AO correction was first conducted on this sample with the intrinsic structures shown in Fig. 5b. With the resolution improvement and signal gain brought by AO, the gap between the exine and body of the grain (the yellow arrow in Fig. 5c) can be clearly distinguished, and porous structures inside the grain can be resolved at higher contrast (Fig. 5d). Next, we acquired the subtracted images (Fig. 5e) of the same lateral section as in Fig. 5a through subtractive imaging with and without AO correction (the wavefront for focus shaping was superimposed with the corrective wavefront in Fig. 5f when AO was applied). To differentiate the features revealed by different image contrasts, the images in Fig. 5e, g, and h are normalized to their respective peak intensities. As shown in Fig. 5g, more pores can be clearly identified in the body of the grain in the subtracted image with AO correction. Although subtractive imaging is proposed to boost the resolving ability of conventional imaging, the quality of the subtracted image can be even worse with the existence of aberrations. As shown in Fig. 5h, without AO correction, more artifacts can be found in the subtracted image than in the conventional image, which makes the blurred boundary indistinguishable from the artifact-corrupted background in the subtracted image. In contrast, while the boundary can be clearly discerned in the conventional image with AO correction, the exine separated from the body of the grain (the cyan arrow in Fig. 5h) can be further resolved at the boundary by subtractive imaging. The results demonstrate the effectiveness of our method for subtractive imaging and suggest that AO correction is necessary to achieve the desired resolution and contrast enhancement.

Aberration correction for structures with temporally varying signals
Instead of retrieving aberrations from the variation in the signal intensity, which is susceptible to signal fluctuations, our approach relies on structure comparison and can be applied to samples with temporally varying signals. To demonstrate this capability, we numerically simulated the AO correction for a sample that consists of two axially separated hollow spheres with identical brightness (Fig. 6a) and applied artificial aberration that is generated in the same way as described in the beads experiment. According to Fig. 6b (No AO), the image of the hollow spheres was severely distorted by the aberration, and neither of the spheres could be discerned. When the aberration was measured with constant signals, isolated spheres were clearly identified and the peak intensity was increased 8-fold after AO correction (Fig. 6b (AO CS)). To simulate signal fluctuations, we multiplied the contributions from different spheres in each acquired image frame by respective factors, which were randomly generated from 0.5 to 2, during each aberration measurement trial. In a representative aberration measurement trial with fluctuating signals, whose modulation sequences consisting of the multiplicative factors for each image to produce the signal fluctuations are shown in Fig. 6c, we observed similar improvement of the image quality as in AO correction with constant signals, while the peak intensity was increased 7.8-fold (Fig. 6b (AO FS)). As shown in Fig. 6d, the features of the hollow spheres recovered in both cases closely matched. In ten trials of AO correction with fluctuating source signals, the peak intensity of the image was increased 7.7-fold on average, which, along with the above results, suggests the robustness of our method to signal fluctuations.
We then applied our method to imaging neurons in the hypothalamus of a living larval zebrafish from transgenic line Ki(th:Gal4-VP16);Tg(UAS:GCaMP6s) with temporally varying fluorescence signals resulting from spontaneous neuronal activity in the brain. As shown in Fig. 7a, the neuronal structures at imaging depths from 262 to 267 μm below the skin were obscured by the aberration. With the neuronal structures axially centered 274 μm below the skin (Fig. 7b), we measured the sample-induced aberration, in which we could identify the presence of astigmatism probably stemming from the cylindricallike surface of the zebrafish larva (Fig. 7c). Correction of the sample-induced aberration led to images with higher resolution and sharper fluorescence signals (Fig. 7d, e) and allowed finer neuronal processes to be better resolved. Next, we recorded the spontaneous activity of the neurons at an imaging depth of 279 μm below the skin. To give a reasonable comparison between images without and with AO correction, image frames without and with AO correction were acquired in an interleaved manner. The maximum intensities of the neuronal structures during calcium imaging without and with AO correction are shown in Fig. 7f. Although bright neurons and fibers could be observed both without and with AO correction, faint structures were only detectable with AO correction. We further monitored the fluorescence signal variation for a few faint features (marked by cyan arrows in Fig. 7f ) during the imaging session and observed more calcium transients after AO correction (Fig. 7g), which suggests the benefit of AO correction for detecting the calcium activity of faint targets.

Conclusions
In summary, we have developed an indirect wavefront sensing method, which uses a virtual imaging scheme and a structural-similarity-based shift measurement to derive aberrations from microscopy images. The presented method could be viewed as a variant of pupil-segmentation-based methods, which have been employed to achieve high-resolution tissue imaging in numerous applications [11-13, 15, 26, 27]. For the conventional approach with single segment illumination [11], the intrinsic structures of biological samples, which are usually not sparsely distributed, present challenges for aberration measurement. In these scenarios, reference beads should be introduced into the samples to assist aberration measurement [12,13], which may not be preferred since this not only complicates the preparation process but also perturbs the native state of the biological samples. Leveraging a virtual imaging scheme, our method is capable of measuring aberrations with more diverse intrinsic structures, and its effectiveness has been demonstrated in both fixed biological tissues and living animals. The full-pupil illumination approach is not limited by the sample structure and can be efficiently implemented via a frequency multiplexing strategy using a high-speed DM [15,26]. However, this method requires the fluorescence signal to be stable during wavefront gradient measurement. Employing a structural-similarity-based shift measurement method, our method can be applied to structures exhibiting fluctuating signals, and the robustness of our method to signal variation has been tested by simulations and aberration measurements with firing neurons. Since our method only requires a commonly used SLM, which is also widely adopted in resolution enhancement techniques, as a wavefront shaper, our method can be directly incorporated with these techniques to image biological tissues with higher resolution and contrast, as shown in AO-incorporated subtractive imaging. Therefore, our method would be a powerful complement for existing pupil-segmentation-based methods that could benefit studies in various fields, such as neurobiology and oncology, by enabling high-resolution imaging within thick tissues.

System configuration
The homebuilt two-photon microscope is depicted in Fig. S2. The two-photon excitation source is a tunable mode-locked titanium sapphire femtosecond laser (Mai Tai DeepSee, Spectra-Physics), with its power controlled by a half-wave plate mounted on a motorized rotation stage (PRM1Z8, Thorlabs) and a polarizing beam splitter cube. The laser (920 nm for zebrafish imaging and 800 nm for the remaining experiments) is expanded to overfill the display panel of a reflective phase-only liquid crystal SLM (X13138-07, Hamamatsu). A second half-wave plate orientates the polarization of the laser for effective wavefront modulation. A pair of lenses (AC254-250-B and AC254-100-B, Thorlabs) relays the modulated light to a scanning unit. The scanning unit is composed of a pair of galvanometers (6-mm aperture, 6215H, Cambridge Technology) conjugated by a pair of lenses (AC508-080-B and AC508-080-B, Thorlabs). A field stop is placed at the intermediate image plane between the SLM and the scanning unit to block undesired diffraction orders. The scanned beam is further relayed by a scan lens (SL50-2P2, Thorlabs) and a tube lens (TTL200MP, Thorlabs) to the back focal plane of a microscope objective (XLPLN25XWMP2, NA 1.05, 25 ×, Olympus). The SLM, galvanometers and back focal plane of the microscope objective are mutually conjugated by the three pairs of lenses. A quarter-wave plate is used to change the polarization of the incident laser into circular polarization. The microscope objective is mounted on an objective scanner (ND72Z2LAQ, Physik Instrumente) for three-dimensional imaging. Along the detection path, two-photon excited fluorescence is separated from the excitation laser and split into two spectral components by dichroic mirrors (FF665-Di02 and FF560-Di01, Semrock) and collected by photomultiplier tubes (H7422-40, Hamamatsu). Bandpass filters (FF02-613/73, FF02-525/40 and FF02-447/60, Semrock, for red, green and blue fluorescence, respectively) are used to purify the collected spectral components. A custom LabVIEW program is used for hardware control and image acquisition. By applying a certain linear phase ramp on the selected subregion of the SLM display, only the light coming from this subregion can pass through the field stop so that the corresponding zonal illumination can be achieved. The computational tasks are implemented with MATLAB scripts and executed on a Dell Precision T7920 workstation (processor: Intel Xeon Gold 5118, graphics card: Nvidia GeForce RTX2080 Ti). The system aberration is corrected before imaging experiments, and images labeled 'No AO' are acquired with system correction.

Wavefront segmentation
We use mask approaches described in [11] to segment the back pupil of the objective and correspondingly the wavefront therein. The pupil is divided into N square subregions (the peripheral nonsquare subregions are not used in our experiments). To measure the gradient of each wavefront segment, each corresponding subregion or mask can be illuminated one at a time, which corresponds to the independent mask approach. Since a single subregion of the pupil may not deliver sufficient signal for aberration measurement, a mask corresponding to contiguous subregions can be simultaneously illuminated, where the gradient being measured can be seen as the average of the gradients of the corresponding constituent wavefront segments. In this case, overlapped masks are used for aberration measurement, which corresponds to the overlapping mask approach. In two of our experiments, we used a 2 × 2 with 1 × 1 stepped overlapping mask approach, in which each mask consists of 2 × 2 subregions and is displaced from its neighboring masks by 1 subregion both horizontally and vertically.

Virtual imaging
According to the layout of the wavefront segments, the point spread functions (PSFs) required for virtual imaging are obtained prior to imaging experiments. For an input vectorial field, ⇀ E i , at the back pupil plane of the objective, the focus field, ⇀ E o , can be obtained by vectorial focus field calculations [28]. For zonal illumination, the amplitude of ⇀ E i is set as 0 beyond the corresponding illuminating aperture. The PSF of two-photon imaging can then be calculated as The volume of the sample structures to be acquired for virtual imaging can be estimated from the excitation volume of each ideal beamlet focus, which we approximate as the volume with at least 1/e 2 of the maximum intensity of the corresponding calculated PSF. Assuming that the effective NA under zonal illumination is NA e , the height of the sample volume can be estimated as H = 1.064 / n − n 2 − NA 2 e ) [29], where λ is the wavelength and n is the refractive index of the immersion medium of the objective. Then, the width of the image volume can be estimated as L = tan (θ)H + l, where θ is the maximum convergence angle and l is the width of the FOV of shifted images, assuming that the FOV is a square. In practice, we acquire a three-dimensional stack of images of the sample structures within the volume axially centered at the central imaging plane, with a width of L ′ and a height of H′, where l < L ′ ≤ L and 0 ≤ H ′ ≤ H. L ′ and H′ are adjusted according to the structure distribution. For example, when the structures are sparsely distributed, H ′ can be determined as 0, which means that only a single image is acquired to construct the reference image and can be seen as the reference image choice in the conventional approach [11]. As the structure distribution becomes denser, larger L ′ and H ′ are required. The imaging parameters used in our experiments are listed in Table S1.
With the acquired image stack as the sample and the calculated PSFs, each reference image can then be constructed from the convolution between the sample and the corresponding PSF. If the sampling intervals of the reference image and the shifted image are different, then linear interpolation is applied to the reference image to match the sampling interval of the shifted image.

Image shift measurement
Given a reference image F and a shifted image G, the goal is to find the shift vector ⇀ s for G such that the similarity between the underlying structures in G and their counterparts in F is the largest with respect to ⇀ s . Considering that the intensity of the same underlying structures could change due to signal variation and deviation of the calculated PSFs from the real ones, the similarity measure should be robust to intensity variation. Inspired by the structural similarity index (SSIM), which is widely used for image quality assessment [30], we adopt the structure comparison term of the SSIM, which is robust to intensity variation, to measure the structural similarity for shift measurement. Since the intensity variation could be spatially variant, local statistics of images are used to calculate the structural similarity, and the structural similarity between image a and image b, both consisting of n × n pixels, is thus given as where a⇀ p and b⇀ p are the patches, both consisting of 11 × 11 pixels, centered at position ⇀ p in a and b, respectively, and where w is a 11 × 11 Gaussian smoothing kernel with a standard deviation of 1.5 pixels that is normalized to unit sum.
The image shift ⇀ s is obtained in two steps. In the first step, both F and G are downsampled k (k = 8) times in both dimensions as F k and G k with the average method. By searching all the patches of the same size as G k in F k for the maximum similarity, a coarse estimate of ⇀ s is obtained as ⇀ s k . In the second step, we generate the corresponding fine estimates ⇀ s ′ from ⇀ s k by For each fine shift ⇀ s ′ , we obtain the corresponding patch of the same size as G in F, namely, F ′ , and the similarity is calculated between F ′ k , which is obtained similarly to F k , and G k . The fine shift that gives the maximum similarity is eventually determined as ⇀ s .

Wavefront reconstruction
Based on the measured gradients of all wavefront segments, a modal reconstruction algorithm and a zonal reconstruction algorithm, which are similar to those described in e , −k/2 ≤ e x , e y < k/2. [22,23], are each used to reconstruct the corrective wavefront. The reconstructed wavefront that gives a higher signal intensity with the same input laser power for the objective is then chosen as the final corrective wavefront. Next, we present the wavefront reconstruction algorithms. Let N denote the number of segments, Σ i denote the region of the i-th segment, i denote the coordinates of the center of Σ i , and ⇀ g i = g x i , g y i denote the measured gradient of the i-th segment.

Modal reconstruction algorithm
In modal reconstruction, the wavefront φ is approximated by the sum of a series of Zernike polynomials such that

Zonal reconstruction algorithm
In zonal reconstruction, the entire wavefront φ is expressed in a segment-wise manner. Let φ i denote the representation of the i-th wavefront segment; then, x, y is the normalized pupil coordinates, ' • ' denotes the scalar product operation, and p i is a constant to be determined. According to [22], we have where E(i) denotes the set of segments that are adjacent to the i-th segment. Then, {p i } can be obtained via an iterative approach or by solving a matrix equation.
When an overlapping mask is used, the wavefront of an overlapping region is set as the average of the covering wavefront segments.

Wavefront correction considerations
The overall aberration measurement time (T m ) consists of two parts: image acquisition (T a ) and aberration calculation (T c ). The image acquisition time T a (N, M) = T 2d (N) + T 3d (M), where T 2d (N) is the time for acquiring N shifted images, and T 3d (M) is the time for acquiring M images in 3D imaging. The aberration calculation time for N wavefront segments is T c (N), which accounts for the construction of reference images, calculation of image shifts and reconstruction of the corrective wavefront. For typical parameters used in our experiments (N = 37, M = 41), the images were acquired at a frame rate of approximately 2 ~ 3 Hz and the computational overhead was less than 0.7 second for each wavefront segment, which results in an overall measurement time with less than 70 seconds. The obtained corrective wavefront is then applied for aberration correction in the formal experiment. As long as the optical aberration is stable in the biological tissues, the corrective wavefront will remain effective. In our experiments, we found that the corrective wavefronts remained effective over the whole experimental duration.

Subtractive imaging
An implementation of subtractive imaging similar to [21] is adopted here. A conventional Gaussian focal spot and a donut-shaped focal spot, which is shaped through 0 ~ 2π vortex phase modulation, are employed to probe the sample. A subtracted image with enhanced resolution and contrast is then constructed by where I c denotes the normalized conventional image acquired with the Gaussian focal spot, I d denotes the normalized image acquired with the donut-shaped focal spot, I s denotes the subtracted image, and γ denotes the weighting factor. Negative values resulting from the subtraction operation are set to zero. In our experiment, γ is set as 0.6, and I c and I d are smoothed by a Gaussian kernel with a standard deviation of 1.5 pixels to suppress noise. When AO is applied, the modulation pattern for desired focus shaping is superimposed with the wavefront correction pattern.

Fluorescent bead preparation
A suspension of 2 μm yellow-green carboxylate-modified microspheres (F8827, Thermo Fisher) was mixed with 1% low-melting-temperature agarose, dispersed by means of a bath sonicator and then mounted onto the bottom well of a glass bottom dish (35 mm dish with a 15 mm bottom well). The bottom well of the dish was fully filled by the sample and sealed with a coverslip.

Brain slice preparation
Mice expressing mRuby3 via virus injection were transcardially perfused with 1X phosphate buffered saline (PBS) followed by 4% paraformaldehyde (PFA). The brains were postfixed with 4% PFA overnight and then washed five times with PBS. Sections of 400 μm thickness were cut on a Leica VT1200S vibrating microtome, and selected sections were mounted onto glass slides and embedded with antifade mounting medium under coverslips.