Laser scanning reflection-matrix microscopy for aberration-free imaging through intact mouse skull

A mouse skull is a barrier for high-resolution optical imaging because its thick and inhomogeneous internal structures induce complex aberrations varying drastically from position to position. Invasive procedures creating either thinned-skull or open-skull windows are often required for the microscopic imaging of brain tissues underneath. Here, we propose a label-free imaging modality termed laser scanning reflection-matrix microscopy for recording the amplitude and phase maps of reflected waves at non-confocal points as well as confocal points. The proposed method enables us to find and computationally correct up to 10,000 angular modes of aberrations varying at every 10 × 10 µm2 patch in the sample plane. We realized reflectance imaging of myelinated axons in vivo underneath an intact mouse skull, with an ideal diffraction-limited spatial resolution of 450 nm. Furthermore, we demonstrated through-skull two-photon fluorescence imaging of neuronal dendrites and their spines by physically correcting the aberrations identified from the reflection matrix.

C onfocal imaging is one of the most widely used configurations in deep-tissue optical microscopy due to its ability to reject signals from out-of-focus planes and other unwanted noise by multiple light scattering 1,2 . The technique has been employed either by itself in confocal reflectance/fluorescence imaging 3 or combined with temporal gating to further attenuate multiple-scattering noise 4 . In the presence of sample-induced aberrations, however, signals containing object information are spread away from the confocal pinhole by the blur of the pointspread function (PSF) [5][6][7][8] . Image contrast is thereby reduced, and resolving power is gradually lost with increasing imaging depth.
The key to resolving sample-induced aberrations and achieving ideal diffraction-limited imaging deep within a scattering medium is to coherently refocus the non-confocal signals, or signals arriving at positions other than the confocal pinhole, back to the confocal detection position. Since multiple-scattering noise arrives at non-confocal positions as well as spatially aberrated signals, it is critical to selectively refocus the aberrated signals, but not the multiple-scattering noise. Up to this point, selective refocus has been extremely difficult in label-free reflectance imaging compared with fluorescence imaging, mainly because the signal and multiple-scattering noise have identical optical frequencies.
The most straightforward approach for focusing non-confocal signals back to the confocal position is to place a wavefront shaping device such as a deformable mirror or liquid-crystal spatial light modulator (SLM) in the illumination and/or detection beam paths. Selective focusing of the aberrated signal is achieved by adjusting the SLM iteratively to maximize the intensity or sharpness of the resulting confocal image [9][10][11][12][13] . These so-called sensorless adaptive optics (AO) methods can identify high-order aberrations, mainly in fluorescence imaging. However, the feedback process is time-intensive because an image acquisition is required for each iteration step. A model-based, modal aberration correction approach might reduce optimization time, but only works on the lowest modes of Zernike polynomials 14 . Furthermore, the iterative optimization can only be initiated when the objects are visible in the first place.
Various approaches have been proposed based on direct wavefront measurements to correct aberrations with the minimal number of image acquisitions; these are referred to as wavefrontsensing AO. The Shark-Hartmann wavefront sensor [15][16][17][18][19] and interference microscopy [20][21][22][23] have been widely used for wavefront recording. In general, sample-induced aberrations in both the illumination and detection beam paths contribute to the measured wavefront. Mathematically, the detection PSF is convolved with the multiplication of the object function and illumination PSF, and this makes it difficult to extract individual PSFs. Previous studies have resolved this double-pass problem in limiting cases when there exist bright, point-like scatterers within the sample, similar to guide stars in astronomy 6 . This changes the multiplication of the object function and illumination PSF to a point source regardless of the distortion of the illumination PSF. However, guide stars are available only in special cases. For example, photoreceptor cells in retina imaging 18 and nonlinear excitation can serve as intrinsic and approximate guide stars 15,16,19,24,25 . Another strategy is to illuminate the planar wave, but not the focused wave, to minimize aberrations in the illumination beam path. Aberrations are then obtained from a reflected wavefront by computationally optimizing image metrics 22,23,26 . The drawback of this approach is its susceptibility to multiple-scattering noise due to the loss of confocal gating. These limitations of existing wavefront-sensing AO approaches are mainly due to incomplete recording of the input-output response of the specimens, which makes the double-pass problem underdetermined.
We hereby propose a laser-scanning reflection-matrix microscopy (LS-RMM) method that records both non-confocal and confocal signals of elastic backscattering from the sample in their phases and amplitudes. Low-coherence interferometry is employed for time-gated detection of the backscattered waves; this excludes multiple-scattering noise arriving at different flight times from the signal waves. Together, these measurements constitute a reflection matrix that quantifies the complete inputoutput response of the light-medium interaction. Conventional optical coherence microscopy (OCM) combining confocal microscopy and optical coherence tomography measures a subset of the reflection matrix. We take advantage of non-confocal signals to unambiguously identify one-way aberrations from highly complex round-trip distortions. The processing step is the transformation of the measured reflection matrix taken for focused illuminations to that for planar illuminations. We then administered a unique algorithm referred to as the closed-loop accumulation of single scattering (CLASS), developed previously to selectively identify and computationally correct sampleinduced aberrations from background multiple-scattering noise [27][28][29] . In particular, we improved the algorithm to correct a large number of correction modes (~10,000 angular modes) at each local area down to 10 × 10 µm 2 . By using this capacity for correcting local high-order aberrations, we demonstrated labelfree, in vivo imaging of myelinated axons underneath an intact mouse skull, an extreme form of aberration medium 30 . Furthermore, by displaying the conjugation of the identified aberration maps on a spatial light modulator in the excitation beam path, we performed two-photon fluorescence (TPF) imaging of the neuronal dendrites and visualized their spines with the spatial resolution of 500 nm, close to the diffraction-limited TPF imaging resolution of 380 nm, over the field of view much larger than the isoplanatic patch.

Results
Experimental schematic of laser-scanning reflection-matrix microscopy. The LS-RMM is built on the same backbone as the OCM system, but with a difference in the detection scheme. A camera was placed on a plane conjugate to the image plane (instead of the confocal pinhole and a photodetector), and a reference wave was introduced to the camera to form off-axis low-coherence interferometry (Fig. 1a, "Methods", and Supplementary Note 1). The time-gated electric-field (E-field) images for waves backscattered from the sample were then retrieved from the interferogram (Fig. 1b-e and Supplementary Note 2). In the absence of aberrations, the obtained E-field image was sharply focused with a diffraction-limited spot size of 450 nm (intensity map, Fig. 1b; phase map, Fig. 1c). On the contrary, the intensity map taken in the presence of an aberrating medium showed a significantly broadened distribution with the peak intensity attenuated by almost two orders of magnitude (Fig. 1d). The corresponding phase map also showed a broadened distribution of meaningful phase values (Fig. 1e).
Image acquisition was conducted using the focused illumination beam scanned in raster mode by the two galvanometer mirrors (GMs) over a given field of illumination (FOI) area on the sample plane with a diffraction-limited scanning interval of λ/ 2α, where λ is the illumination wavelength, and α is the numerical aperture (NA) of the objective lens. The FOI is equivalent to the field of view of the conventional OCM. For each illumination point r i at the sample plane, the E-field image E cam (r cam ; r i ) was recorded in the de-scanned camera coordinate frame with the position vector r cam . In contrast to previous wavefront-sensing AO, where a wavefront detector was located at a conjugate pupil plane 20,31 , the present scheme offers the best spatial distinction between the signal and multiple-scattering noise 32 . Here, the detection area at the camera, called the field of detection (FOD), was adjusted to be wide enough to capture the entire profile of the reflected wave broadened by the sample-induced aberrations. The number of correction modes N c is determined by the number of free modes within the FOD. As discussed in detail in "Methods", in a typical FOD of 50 × 50 µm 2 , we can correct up to 10,000 angular modes of aberrations in the pupil. The choice of optimal sizes for the FOD and FOI is made in consideration of PSF broadening, the intensity ratio between single-and multiplescattered waves, and image acquisition time. The wider the FOD used, the higher order of aberrations can be measured and corrected, with a concomitant increase in data acquisition time (Supplementary Note 5), and the more multiple-scattering noise is captured. We demonstrated the hardware correction of the identified aberration by inserting an SLM (X10568-02, Hamamatsu) at the plane conjugate to the objective lens in the illumination beam path. The SLM serves as a flat mirror for reflection matrix recording and image acquisition of label-free reflectance imaging. The device was used for multi-photon imaging to physically compensate for the sample-induced aberrations identified by the CLASS algorithm.
The LS-RMM proposed here is different from previous reflection-matrix-based approaches in its illumination and detection configurations 27,28,33 . Our earlier studies used planar wave illuminations, while the present LS-RMM uses a focused illumination. The main benefit therein is the enhanced signal to multiple-scattering noise ratio. In planar wave illumination, multiple-scattering noise generated by illumination over a wide area is summed at each detection pixel. On the contrary, signals generated by a focused illumination compete only with the multiple-scattering noise generated by the given focused illumination itself 32 . Other approaches employing focused illumination usually put a detector in the pupil plane where signals are spread out over the detection area 20,31,34 . In contrast, we placed a wavefront detector at the image plane where signals are more likely to be focused relative to multiple-scattering noise, providing an optimal spatial distinction between the signal and multiple-scattering noise. This technique ensures the optimal use of the detector's dynamic range, resulting in enhanced sensitivity to multiple-scattering noise (Supplementary Note 4). Furthermore, the matrix acquisition time can be shortened in the new configuration. The view field for recording non-confocal signals needs to be wide enough to capture the entire PSF broadened by the sample-induced aberrations. We can adaptively adjust this view field at the camera with the increase of imaging depth in such a way to optimize matrix acquisition speed. This is especially critical for in vivo imaging where the motion of the living specimen can undermine the image reconstruction process.  Another important benefit of focused illumination geometry lies in its direct compatibility with existing imaging modalities such as OCM and multi-photon microscopy (MPM). The AO system can be integrated with these conventional imaging modalities by inserting wide-field, low-coherence interferometry in the detection plane, allowing the combination of hardware-adaptive optics (HAO) and computational adaptive optics (CAO) 21 in complement to one another. We realized HAO by displaying the conjugation of the aberration maps identified from the reflection matrix by simply inserting the SLM in the illumination beam path and demonstrated AO multi-photon imaging through an intact mouse skull. We would like to emphasize that the realization of LS-RMM in a focused illumination configuration has been technically challenging due to shot-to-shot phase fluctuations during scanning, especially for in vivo imaging applications. We maintained phase-referencing between different illuminations by improving the mechanical stability of the interferometry and laser-scanning system.
Image processing and aberration correction. In the case of optical coherence reflectance imaging, sample-induced aberrations of the incident wave on the way to the sample and those of the reflected wave on its way back to the detector must be separately identified. This necessitates the recording of a reflection matrix and the application of the CLASS algorithm. In this section, we described how the sequence of images taken by LS-RMM was processed to construct a time-gated reflection matrix in the space domain. For a given focused illumination at a position r i in the sample plane, the time-gated E-field of the reflected wave E cam (r cam ; r i ) recorded at the position r cam in the camera plane is as follows (see Supplementary Note 3 for the coordinate system): Here, P i (r; r i ) is the illumination E-field PSF at a position r in the sample plane, which is broadened with respect to r i due to sample-induced aberrations. O(r) is the object function given by the amplitude reflectance of the target object. P o (r cam ; r i ) describes the detection E-field PSF of the returning waves from the sample plane to the camera plane. E M (r cam ; r i ) represents a speckle field from the time-gated multiple-scattered waves by the scattering layer located above the sample plane. Multiplescattered waves whose flight times are either longer or shorter than that of the signal wave are rejected by the time-gated detection.
We next needed to obtain the E-field map of the reflected wave E lab (r o ; r i ) in the laboratory frame coordinate r o to identify the sample-induced aberrations. Since E cam (r cam ; r i ) is acquired after the de-scanning action by the GMs, r cam and r o have the relation For demonstrating the working principle of the LS-RMM system, we performed reflectance imaging of a custom-made Siemens star target placed under a 600-µm-thick, rough-surfaced plastic layer exhibiting strong aberrations (Fig. 2a). A set of E-field images E lab (r o ; r i ) for all illumination spots r i within an FOI of 40 × 40 µm 2 was recorded (Fig. 2b). A time-gated reflection matrix R(r 0 ; r i ) in the space domain was constructed by assigning E lab (r o ; r i ) to each element of R in such a way that the r i and r o are the column and row indices of R, respectively (Fig. 2c). The acquired reflection matrix R is represented as Here, O is a diagonal matrix containing the target's object function in the diagonal elements; O(r; r) = O(r). P i is the illumination PSF matrix whose element is given by P i (r; r i ). The detection PSF matrix P o consists of P o (r o ; r), which is theoretically identical to the transpose of P i in epi-detection geometry but can differ in the presence of a slight mismatch between illumination and detection optics. M is the multiplescattering matrix composed of E M (r o ; r i ). The conventional OCM image can be obtained from the main diagonal elements of R as they correspond to the confocal signals (Fig. 2d). The OCM image was blurred and severely distorted due to pronounced sample-induced aberrations contained in P i and P o . We then identified and removed P i and P o in R to reconstruct an aberration-free object image.
To apply the CLASS algorithm, we converted the position basis of the reflection matrix R into the spatial frequency basis (k space) by applying a Fourier transform operator F and the inverse operator F −1 to the output and input sides of R, respectively:R ¼ FRF À1 . The reflection matrixR on the spatial frequency basis (Fig. 2e) can be expressed as Here,Õ ¼ FOF À1 is a circulant matrix whose elements are given by the spatial frequency spectrum of the object functioñ Þwhere, k i and k o are the input and output transverse wavevectors, respectively.P i ¼ FP i F À1 and P o ¼ FP o F À1 are the transmission matrices for a plane wave propagating through the illumination and detection pathways, respectively. For the area where the PSF is shift-invariant, which is the case within an isoplanatic patch, the relation P i (r; r i ) = P i (r − r i ) is valid. Then, P i takes the form of a Toeplitz matrix, which setsP i as a diagonal matrix whose diagonal element consists of exp[iϕ i (k i )]. Likewise,P o is a diagonal matrix with the diagonal element given by exp[iϕ o (k o )]. Here, ϕ i (k i ) and ϕ o (k o ) are angle-dependent phase retardations induced by the scattering layer for the illumination and detection pathways, respectively.
We then applied the CLASS algorithm 27,28 toR, which separately identifies ϕ i (k i ) and ϕ o (k o ) such that the total intensity of the confocal image, reconstructed from an aberrationcorrected reflection matrix, is maximized (see "Methods"). The conjugations of the pupil functions were applied to the original matrixR to obtain the aberration-corrected matrix,R c ¼P * oRP * i , which was then transformed back to the position-space matrix R c (Fig. 2f). The retrieved ϕ i (k i ) and ϕ o (k o ) were almost identical due to optical reciprocity (Fig. 2g, h). The number of angular modes N c used for aberration correction in the pupil was about 6200, which is set by the FOD size of 40 × 40 µm 2 . The magnitudes of the diagonal elements of R c were greatly enhanced owing to the compensation of the sample-induced aberrations, and those of the off-diagonal elements were reduced in comparison to the original matrix R (Fig. 2c). Compared to the original OCM image (Fig. 2d), the intensity of the aberration-corrected image (Fig. 2i) acquired from the diagonal elements of R c was increased by an average of about 30 times. The spatial resolution of the corrected image was estimated to be 450 nm, identical to the ideal diffraction limit.
We would like to emphasize that a new algorithm was added to the original CLASS algorithm in imaging through a mouse skull to correct local position-dependent aberrations ("Methods" and Supplementary Note 6). In complex scattering layer, the assumption that P i r; r i ð Þ ¼ P i r À r i ð Þ is valid only within the small isoplanatic patch. This means that aberrations cannot be corrected by a single pair of ϕ i (k i ) and ϕ o (k o ). We developed an algorithm to identify multiple ϕ i (k i ) and ϕ o (k o ) functions for individual isoplanatic patches. The algorithm is special in that the number of correction modes is maintained at the number of angular modes in full FOD (~10,000) even when the image analysis area is reduced to the isoplanatic patch size of 10 × 10 µm 2 . This is contrary to conventional computational AO, where the number of correction modes is reduced in proportion to the image analysis area.
Physical aberration correction by using SLM. The LS-RMM system was designed not only to identify high-order aberrations induced by the medium but to physically compensate for them to generate a near-diffraction-limited focus on the sample. To this end, the SLM in the experimental setup (Fig. 1a) was used to display the phase conjugation of the identified input aberration map (Fig. 2g)  Each image was converted to a column vector and assigned to its corresponding column in R. d OCM intensity image constructed from the main diagonal of R in c before aberration correction. Scale bar, 10 µm. e Reflection matrixR in spatial frequency space. f Aberration-corrected reflection matrix R c converted from R after application of CLASS algorithm. g, h Phase maps ϕ i (k i ) and ϕ o (k o ) for aberrations in illumination and detection pupils retrieved by the CLASS algorithm, respectively. The radii of the maps correspond to a numerical aperture of 1.0. The number of modes N c used for aberration correction in the pupil was about 6200. i Aberration-corrected CLASS intensity image obtained from the main diagonal of R c in f. j Schematic of hardware wavefront correction. Conjugate of the illumination pupil phase map in g displayed on the SLM in Fig. 1a to physically compensate for the aberrations. k OCM intensity image after hardware correction (HC) of aberrations. l Intensity image of reflected PSF measured at the camera after physical aberration correction by SLM (HC PSF). Scale bar, 10 µm. m Line profiles of the PSFs obtained without wavefront correction (black), after computational wavefront correction by the CLASS algorithm (blue), and after hardware wavefront correction by SLM (red). NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19550-x ARTICLE NATURE COMMUNICATIONS | (2020) 11:5721 | https://doi.org/10.1038/s41467-020-19550-x | www.nature.com/naturecommunications waves incident to and outgoing from the sample all together (Fig. 2j). With this hardware correction in place, we could obtain a clean OCM image (Fig. 2k) with no need for computational correction by the CLASS algorithm. The image quality was comparable to that seen in the CLASS image (Fig. 2i), which confirms that a sharp focus was physically generated at the plane of the target object. The PSF came into sharp focus, with a significantly increased peak intensity, after hardware wavefront correction by the SLM (Fig. 2l). We assessed the quality of the aberration correction by measuring the intensity profile across the center of the PSF (red curve, Fig. 2m). The original PSF of the reflected wave (black curve, Fig. 2m) was highly speckled with increased width. The enhancement by the Strehl ratio, measured by the ratio of the peak intensities after and before aberration correction, was 26, and the width of the focus approached the diffraction-limited spatial resolution of 450 nm. Peak height was about 50% of the computational aberration correction by the CLASS algorithm (blue curve, Fig. 2m), due to SLM limitations in shaping steeply varying aberrations, especially in the high spatial frequency range. The width of PSF, however, was almost identical to that of the computational correction.
In vivo imaging of a mouse brain through an intact skull. We tested the capacity for the CLASS algorithm to correct extreme, high-order aberrations, by attempting to image neuronal structures in the mouse brain through an intact skull. The mouse skull is composed of fine microstructures that induce severe optical aberrations as well as strong multiple-scattering noise. Thus far, only AO TPF fluorescence imaging 10,24,35 or three-photon excitation imaging 36 with excitation wavelengths greater than 1300 nm have been able to visualize mouse brain structures through the skull; label-free reflectance imaging has not yet been able to do so. This is because the round-trip aberrations jointly deteriorate the image quality, while input aberration matters most in multi-photon microscopy. The mouse skull consists of several layers of microstructures; these reduce the size of the isoplanatic patch down to the degree that confounds conventional AO methods 30 . In this study, we made significant improvements to the CLASS algorithm to correct aberrations locally with an isoplanatic patch size as small as 10 × 10 µm 2 while maintaining the number of correction modes to~10,000 modes (Supplementary Note 6), ameliorating the extreme aberrations caused by the skull. Figure 3 shows CLASS images of myelinated axons in the first cortex layer after different types of sample preparation. We first conducted imaging through a thinned skull (Fig. 3a-d), wherein compact bone and spongy bone were removed from the skull of a 7-week-old mouse to reduce its thickness to approximately 40 µm (Fig. 3a). We recorded a time-gated reflection matrix R with a FOD of 30 × 30 µm 2 for the FOI of 150 × 150 µm 2 at a depth of 72 µm from the upper surface of the skull. In a conventional OCM image (Fig. 3b) reconstructed by the diagonal elements of R, the shapes of microvessels (blue arrowheads) and the myelinated axonal fibers (yellow arrowheads) appeared blurred. The effective diameter of a confocal pinhole used in the OCM image was set to 1 Airy unit. To correct anisoplanatic aberrations over a large field of view, we divided the entire view field into 18 × 18 regions and applied the CLASS algorithm to the individual subregions. Both microvessels and myelinated axons were visible with improved sharpness and signal strength in the aberration-corrected image (Fig. 3c). The retrieved aberrations varied spatially, but the degree of aberration was mild (Fig. 3d) due to the thinning of the skull. The size of each subregion was approximately 11 × 11 µm 2 , including margins that overlapped with adjacent subregions. The number of correction modes for each aberration map, set by the FOD, was about 3500.
We then conducted imaging through an unaltered, intact skull excised from an 8-week-old mouse. In the 3D image reconstructed by SHG signals detected at the PMT, compact bone, spongy bone, and meninge could be discretely identified (Fig. 3e). Using the SHG image, the thickness of the skull was measured to be 90-110 µm. The time-gated reflection matrix R was measured with a FOD of 50 × 50 µm 2 for the FOI of 150 × 150 µm 2 at a depth of 150 µm below the upper surface of the skull (red dotted box in Fig. 3e). OCM (Fig. 3f) and CLASS images (Fig. 3g) were obtained along with local aberrations (Fig. 3h). Compared to the thinned skull, the myelinated axons were virtually invisible in the OCM image because the thicker skull induced more complex aberrations. Our advanced CLASS algorithm was able to compensate for these increased aberrations and reconstruct myelinated axons at a high spatial resolution and contrast (Fig. 3g). The smallest thickness of the myelinated axons was measured to be 450 nm, confirming that the ideal diffractionlimited resolution was recovered. The aberration maps identified by our advanced CLASS algorithm with the same patch size as the thinned skull show increased complexity, and the correlation between aberration maps of adjacent regions was reduced. The number of correction modes in each map was about 10,000. From the obtained pupil aberration maps, it is estimated that the PSF width of the aberrated single-scattered waves is about 6-8 µm in full width at half maximum, which causes a reduction in the peak intensity of the single-scattering signal in the confocal spots by a factor of~400. By comparing the PSFs before and after aberration correction, we estimate that the ratio of single-scattering signal to time-gated multiple-scattering background noise at confocal points was initially about 0.08, much smaller than 1, before the aberration correction. This explains why the conventional OCM failed to achieve high-resolution imaging of the mouse brain. The CLASS algorithm selectively refocused the aberrated singlescattering signals back to the confocal points to raise the singlescattering intensity by a factor of~400. This made the singlescattering intensity larger than time-gated multiple-scattering noise by about 30 times after the aberration correction and enabled us to identify individual myelinated axons with diffraction-limited resolution (see Supplementary Note 7 for detailed PSF analysis).
Finally, we conducted in vivo imaging of myelinated axons through the intact skull of a living mouse (Fig. 3i-m). After removing the scalp from an 8-week-old mouse (Fig. 3i), a 100µm-thick cover glass was glued on the surface of the intact mouse skull. A skull holder was attached to the mouse head with dental cement, which was then firmly fixed onto the custom-built stage. The thickness of the skull was measured to be~120 µm. We measured the time-gated reflection matrix R at a depth of 200 µm below the skull with the FOD of 30 × 30 µm 2 and obtained an OCM image from the diagonal elements (Fig. 3j). The aberrations were pronounced, such that no structures could be resolved in the OCM image. We then applied the CLASS algorithm to a 30 × 30 µm 2 sized area (blue dotted box, Fig. 3j) and obtained an aberration-corrected image (Fig. 3k). The reconstructed image shows signatures suggesting the existence of certain structures, but they could not be well resolved because the size of the FOI was larger than the isoplanatic patch size of the skull.
We addressed the local aberrations by dividing the area into 2 × 2 subregions and applied the CLASS algorithm to each subregion (Fig. 3l). Bright fiber structures began to emerge in the reconstructed image, and the retrieved aberration map in Fig. 3l also showed more complex phase structures. We further reduced the size of the subregion to 10 × 10 µm 2 and reconstructed the CLASS image (Fig. 3m). Myelinated axons were resolved without distortion following the aberration correction. Aberrations were complex, and the identified aberration maps resembled the NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19550-x ARTICLE speckle pattern. Of note, the application of the CLASS algorithm to a subregion size greater than 10 × 10 µm 2 failed to resolve the fine structures of myelinated axonal fibers, suggesting that the isoplanatic size of the skull was about 10 × 10 µm 2 . When the subregion is larger than the isoplanatic patch, local aberrations are averaged in such a way that only the slowly varying aberrations in the pupil plane can be corrected; as such, the aberration map in Fig. 3k is much smoother than that in Fig. 3m. The isoplanatic patch size can systematically be found by finding the size of the subregion at which the intensity of the aberrationcorrected image is maximum (see Supplementary Note 8). To the best of our knowledge, these results represent the first experimental demonstration of label-free reflectance imaging through an intact mouse skull at an ideal diffraction-limited resolution.
Near-diffraction-limited TPF imaging through an intact mouse skull. The LS-RMM can directly be combined with multiphoton microscopies such as TPF microscopy and secondharmonic generation (SHG) microscopy to recover their diffraction-limited spatial resolution in deep-tissue imaging. The experimental setup needs little modification, and all needed is to place a PMT right behind the dichroic mirror for detecting the nonlinear emissions from the specimens (Fig. 1). Here, we demonstrated near-diffraction-limited TPF imaging of neuron's dendrites through an intact skull. We prepared a fixed wholemouse brain with an intact skull from a 4-week-old transgenic mouse that expresses enhanced green fluorescent proteins (EGFP) at the neuronal membranes (see "Methods" for sample preparation). We first performed CLASS imaging to obtain the aberrations induced by the mouse skull. The objective focus was set to 125 μm below the upper surface of the skull, whose thickness was about 85 μm. Conventional OCM failed to visualize any myelin due to complex aberrations by the skull (Fig. 4a). We divided the entire field of view into 14 × 14 subregions and applied the CLASS algorithm for each subregion to recover fine myelination structures therein (Fig. 4b) and to obtain the aberration map at the corresponding subregion (Fig. 4c). The number of correction modes in each aberration map was 9880. Similar to the demonstration in Fig. 2j-m, we were able to physically correct the skull-induced aberrations by displaying the phase conjugations of the aberration maps on the SLM. Since the aberration varied depending on the position, we made a physical correction for each subregion at a time. As a point of reference, a conventional TPF image was taken without aberration correction. Figure 4d shows a maximum intensity projection (MIP) of a volume image obtained in a depth range of 119-135 μm with 0.5μm axial spacing over the same field of view like that in Fig. 4a. The contours of neuronal dendrites are blurry, and their microstructures are invisible due to the significant broadening and attenuation of the excitation PSF. For the subregion marked by the yellow dotted box in Fig. 4d, we applied the hardware correction using the aberration map obtained from the same subregion in Fig. 4b and conducted TPF imaging. The TPF image with hardware aberration correction (Fig. 4e) showed 19 times increased in fluorescence intensity and enabled the recovery of the fine dendritic spines that were invisible before the correction. Notably, only the structures in and around the subregion were corrected properly, supporting that the isoplanatic patch size is about the same as the subregion indicated by the yellow box, which is 10 × 10 μm 2 . For the full mapping of the neuronal dendrites, we conducted TPF imaging for each subregion indicated by the dotted gray box in Fig. 4b by displaying the phase conjugation of the associated aberration map in Fig. 4c. Multi-depth TPF images were obtained by scanning the imaging depth over the same depth range in Fig. 4d, where the imaging depth was scanned by adding a defocus phase map to the SLM. Figure 4f, g show MIP images for a depth range of 113 ± 1.5 μm before and after aberration correction, respectively. Figure 4h, i are the same as Fig. 4f, g, respectively, but for a depth range of 122 ± 1.5 μm. The corrected view field was expanded as expected, and dendrites and spines at two different depths were resolved. The measured width of the dendritic spine was as small as 500 nm, close to the diffraction limit.
TPF imaging that we demonstrated here can be considered a type of wavefront-sensing AO. LS-RMM serves as a tool to measure wavefront distortion by scattering tissues, similar to Shack-Hartmann wavefront sensors and other interferometric microscopy techniques. The uniqueness of our approach is that LS-RMM can identify extremely severe skull-induced aberrations without using either external guide stars or nonlinear fluorescence excitation. This benefit comes from the recording of the time-resolved reflection matrix and the use of the CLASS algorithm extracting the one-way aberration. In LS-RMM, the identification of wavefront aberrations is based on the intrinsic reflectance contrast. Myelinated axons are a good source of intrinsic contrast as they are spread throughout the mouse brain. It is not always necessary for the myelinated axons to be around the region of interest for multi-photon imaging. Even for the area where there are no axons such as the subregion indicated by the white arrowhead in Fig. 4b, we were able to identify the aberration map and obtain near-diffraction-limited TPF images. This is because noise-like speckle patterns backscattered from irregular tissue structures at the focal plane are also singlescattered waves contributing to the CLASS algorithm.
Adaptive optics by LS-RMM is beneficial over existing AO modalities that rely upon feedback optimization of multi-photon fluorescence image in that photobleaching effect is negligible during the measurements of the aberrations as it does not require high excitation power to record the reflection matrix. We used an excitation power of 21 mW at the sample and a pixel exposure time of 100 μs for TPF imaging, whereas an excitation power of 3 mW was used for recording the reflection matrix. Another important benefit is that recording of one reflection matrix over the area of 150 × 150 μm 2 enabled us to obtain a 15 × 15 number of aberration maps for as many different isoplanatic patches. In the iterative feedback AO, multiple iterations of optimization need to be conducted for each isoplanatic patch.

Discussion
We developed a laser-scanning reflection-matrix microscopy (LS-RMM) and streamlined the closed-loop accumulation of singlescattering (CLASS) algorithm for correcting local high-order aberrations even in the presence of strong multiple-scattering noise. The system is highly compatible with the conventional laser-scanning confocal microscopy as it shares the same backbone, but it is distinct from conventional confocal imaging in that backscattering signals arriving at both non-confocal and confocal positions are detected. We constructed a time-gated reflection matrix in the space domain and corrected sample-induced wavefront aberrations without requiring guide stars. In particular, we improved the aberration correction algorithm to correct up to 10,000 angular modes of aberrations locally varying at every 10-µm interval in the sample plane. We thus demonstrated labelfree reflectance imaging of cortical myelination with extremely complex local aberrations caused by the intact skull of a living mouse. Since the aberration correction takes place computationally in the coherent reflectance imaging, the recording of the reflection matrix itself is the end of the imaging session. This is particularly beneficial for in vivo imaging because we have ample time to correct severe position-dependent aberrations after the data acquisition. Still, the data acquisition speed is critical for in vivo imaging because even small perturbation for the reflection matrix can undermine the image reconstruction, especially in the case of severe aberrations. Reflectance imaging is more challenging than fluorescence imaging in that the contrast of the reflectance imaging of living tissues is intrinsically low. Furthermore, the round-trip aberrations jointly deteriorate the point-spread-function while only the aberration in the excitation beam matters most in the fluorescence imaging. Our LS-RMM addressed all these difficulties for us to realize the first throughskull reflectance imaging to date.
The working depth of current implementation is set by the weak reflectance contrast of myelinated axons relative to the multiple-scattering noise. LS-RMM can employ any excitation wavelength, similar to confocal reflectance imaging. Therefore, the use of longer excitation wavelengths can help to increase imaging depth as the multiple-scattering noise is to be reduced 37 .
Alternatively, the use of visible wavelengths will lead to increasing the image contrast and spatial resolution, which will open the possibility of investigating detailed myelin pathologies through intact skull 38 . Regarding the aberration correction algorithm, previous reflection-matrix approaches often use singular value decomposition to cope with multiple-scattering noise and/or aberrations 31,34,39,40 . In the context of optical imaging, it was demonstrated to the extent that a sharp image was recovered for a few-micron-sized highly reflecting resolution target hidden under a scattering and aberrating tissue 34 . On the contrary, our proposed algorithm works faithfully for the recovery of diffractionlimited resolution of weakly reflecting biological specimens embedded within biological tissues under a rat skull.
In addition to the reflectance imaging, we physically corrected the skull-induced aberrations by using an SLM in the excitation beam path and generated a sharp optical focus within brain tissues under the skull. This physical aberration correction was applied to the through-skull multi-photon imaging. We demonstrated TPF imaging of neuronal dendrites and their minute spines with the spatial resolution of 500 nm, close to the diffraction-limited resolution of 380 nm. Notably, multiple aberration maps were displayed in sequence for individual 10 × 10 μm 2 subregions to recover dendrites over the field of view much larger than the isoplanatic patch size set by the skull. Furthermore, we also demonstrated the near-diffraction-limited SHG imaging of collagen fibers underneath an intact mouse skull (Supplementary Note 9). In LS-RMM, the identification of wavefront aberrations is based on the intrinsic reflectance contrast of targets. As such, it does not require fluorescent labeling and high excitation power, contrary to existing AO modalities that rely upon multi-photon fluorescence feedback signals. Another critical advantage of LS-RMM is that the computational process to retrieve both the wavefront aberrations and the aberration-corrected reflectance image takes place after the recording of the reflection matrix. The hardware aberration correction capability of LS-RMM can also be added to other imaging modalities such as super-resolution imaging 41,42 and coherent Raman imaging 43,44 to enhance their imaging depth in deep-tissue imaging. The practical limitation of our approach for in vivo multi-photon imaging is the computation time for the CLASS algorithm. It does not take long to process the algorithm for mild aberrations from typical brain tissues. However, it is time-consuming to process the skull-induced aberrations with an extremely small isoplanatic patch, which could be improved through the use of graphics processing units (GPUs) or clusterbased hardware; algorithm improvements in the future could also shorten this step.

Methods
Experimental setup in Fig. 1a. A mode-locked Ti:Sapphire laser (center wavelength, 900 nm; bandwidth, 25 nm; repetition rate, 80 MHz) was used as a shortcoherence-length light source (coherence length, 30 µm). The laser beam was split into sample and reference beams by a beam splitter (BS1), and they were recombined by another beam splitter (BS2) to form Mach-Zehnder interferometry. The sample beam was relayed via two GMs and focused on the sample plane by a water immersion objective lens (OL, Nikon, ×60, NA 1.0). The GMs were used to rasterscan the focused illumination beam on the sample as in a conventional confocal microscope. The backscattered wave from the sample was collected by the same OL and traveled through a reciprocal path. The wave was then de-scanned by the GMs and delivered to a camera (s-CMOS, pco.edge 4.2, PCO AG) placed at the conjugate image plane. The reference beam was sent through multiple pairs of relay lenses, an optical delay line (not shown), and a diffraction grating (DG). The firstorder diffraction of the reference beam by the DG was selected for by an iris diaphragm (I) and delivered to the camera as a plane wave. Consequently, its phase front was tilted with respect to the camera plane while its temporal pulse front remained parallel to the camera plane. The interference between the sample and reference waves formed an interferogram within the temporal gating width corresponding to 15 µm in length, half the coherence length of the light source.
The acquisition time for the full reflection matrix. The frame rate f cam of the s-CMOS camera used here is given by , wherein the experiment, the coefficient r was 40,000 [1/µm]. The number of correction modes N c in a circular pupil with a numerical aperture, α is given by is the diffraction-limited resolution. For example, an FOD size of 50 × 50 µm 2 covers approximately 10,000 correction modes for α = 1 and λ = 900 nm. The number of images to be recorded for complete measurement of a reflection matrix for a given FOI is derived by N s ¼ FOI=δ 2 d . The total acquisition time for a full reflection matrix is then determined by T ¼ N s =f cam . For the FOI of 50 × 50 µm 2 and FOD of 50 × 50 µm 2 , it takes about 15 s to measure the reflection matrix. At a smaller FOD size of 16 × 16 µm 2 containing about 1000 correction modes, the data acquisition time is reduced to around 5 s. The use of a high-speed camera along with optimization of the algorithm could potentially reduce this acquisition time to below 1 s.
Aberration correction using the CLASS algorithm. The CLASS algorithm is an iterative method for finding the aberration maps ϕ i (k i ) and ϕ o (k o ) independently fromR k o ; k i ð Þ. In nth iterative step, the trial solutions for the aberration maps, Here, Δk = k ok i is the momentum difference, andR n ð Þ Δk; k i ð Þand Advanced algorithm for correcting spatially varying high-order aberrations.
To correct local aberrations, it is necessary to analyze the FOI equal to or smaller than the isoplanatic patch given by the scattering medium. In conventional computational AO and the previous CLASS algorithm, dividing the entire FOI into small subregions accompanies the simultaneous reduction of FOI and FOD in each subregion. This results in the reduced N c , making it difficult to correct complex aberrations in small isoplanatic patches. In this study, we developed an algorithm that can find local aberrations while N c is maintained to be large enough. In LS-RMM, we can set FOD independently of FOI to make use of the information outside the subregion. For example, in our in vivo imaging through the intact mouse skull (Fig. 3m), we set the FOD (30 × 30 µm 2 ) larger than FOI (10 × 10 µm 2 ). Consequently, the reflection matrix R sub (r o ; r i ) for each subregion became an N o × N i rectangular matrix, where N o is set by the number of pixels in ffiffiffiffiffiffiffiffiffiffi FOD p þ ffiffiffiffiffiffiffiffi FOI p À Á 2 and N i is set by the number of pixels in FOI. Therefore, this local reflection matrix contains only 390 input modes set by FOI, but it contains N c = 3500 output modes set by the FOD. In our advanced CLASS algorithm, we added zero columns to R sub (r o ; r i ) to convert it into an N o × N o square matrix. We then transformed it intoR sub k o ; k i ð Þhaving the same grid spacing δk ¼ 2π= ffiffiffiffiffiffiffiffiffiffi FOD p þ ffiffiffiffiffiffiffiffi FOI p À Á for both the input and output pupil planes. In the nth iteration, the output aberration map ϕ ðnÞ o k o ð Þ was identified from Eq. (5) for all the N c modes supported by the FOD. Based on the optical reciprocity, we applied ϕ ðnÞ o Àk o ð Þ to ϕ ðnÞ i k i ð Þ to findR n ð Þ sub k o ; k i ð Þin Eq. (6). By repeating this process, we were able to correct both the input and output aberrations for all the N c angular modes even though the FOI is reduced to the size as small as an isoplanatic patch. Due to this capacity, we could identify locally varying high-order aberrations (Fig. 3m). For detailed analysis, see Supplementary Note 6.
Sample preparation for ex vivo imaging. All animal experiments were approved by the Korea University Institutional Animal Care & Use Committee (KUIACUC-2019-0024). Three-to eight-week-old C57BL/6 mice or Thy1-EGFP line M (Jackson Labs #007788) mice were deeply anesthetized with an intraperitoneal injection of ketamine and xylazine at 100 and 10 mg/kg, respectively, and decapitated. The scalps were completely removed, and the heads were fixed with 4% paraformaldehyde at 4°C for 1−2 days. Heads were then washed with phosphate-buffered saline (PBS) in triplicate. For imaging through the intact skull, an unaltered head was affixed to a 35-mm dish using cyanoacrylate and filled with PBS. For imaging through the thinned skull window, compact and spongy bones were ground from a fixed head. The 3-4-millimeter area in the center of the parietal bone was thinned using a 0.5-mm surgical carbide bur in a high-speed dental drill. Afterward, the bone was smoothed with a stone bur and a polishing bur to a thickness of 30−40 μm.
Preparation of intact skull window for in vivo imaging. We followed the procedure outlined by Wang et al. 36 with slight modifications. Briefly, mice (6-8 weeks old, 22−28 g) were temperature-controlled and anesthetized with isoflurane (1.5 −2% in oxygen to maintain a breathing frequency of around 1 Hz). Dexamethasone (1 mg/kg) was administrated intramuscularly on the 2 consecutive days after the surgery to minimize swelling at the surgery site. Eyes were covered with eye ointment during surgery and imaging. The hair was removed with scissors and Nair hair remover. The scalp was removed to expose the bregma, lambda, and parietal plates. Sterile saline was applied to the skull, and connective tissue remaining on the skull was gently removed. A sterile, round, 4 mm (diameter) coverslip (#1, Warner Instruments, USA) was attached to the center of the parietal bone using ultraviolet-curable glue (Loctite 4305). A custom-made metal plate was attached to the skull with cyanoacrylate for head fixation, and the exposed part of the skull was covered with dental cement (Dentsply DeTrey GmbH, Germany).
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data that support the findings of this study are available from the authors on reasonable request.

Code availability
The code used in this study is available from the corresponding author upon reasonable request.