Noninvasive focusing of light and imaging of objects embedded in scattering tissue are crucial for both basic biological research and clinical applications. Tissue scattering is the primary limitation for deep tissue focusing and imaging. Advanced microscopy techniques, such as confocal microscopy and multiphoton microscopy, enable noninvasive deep tissue imaging. However, as imaging depth increases, the number of ballistic photons available for excitation at a target location decreases exponentially. Eventually, the background intensity becomes comparable to the signal in non-sparsely labeled samples [1].

To increase fluorescence signal intensity in two-photon microscopy, adaptive optics has been employed to correct low-order aberrations [2,3]. Recent cutting-edge research [4,5] demonstrates that combining adaptive optics with multiphoton microscopy can enhance the fluorescence signal by several folds. A spatial light modulator (SLM) or a deformable mirror, positioned at the conjugate plane of the objective lens’s back aperture, is used to modulate the phase of light in the Fourier domain. However, the patterning speed of SLM is relatively slow (<1kHz), limiting its application to correct dynamic scattering [6].

Numerous techniques have been developed to correct highly spatially varying scattering [7–11]. However, most of these techniques are designed for transmission mode, which limits their applications for imaging objects embedded within the scattering tissues of living animals. Among the reflection mode techniques, many are post-processing methods [12–16], including our previous work [17–21], which require significant time for image reconstruction ranging from several minutes to over an hour. Unlike post-processing techniques that correct scattering in images offline, several techniques achieve active correction of scattering by wavefront shaping using an SLM in reflection-mode multiphoton microscopy systems [22–25]. The wavefront shaping is based on interferometry measurements. Excitation light is split into two beams, either in a common path [22,24] or in a divergent path [23,25]. One beam is steady, while the other beam is phase-shifted using either an SLM or a galvo-scanner. The interference of these two beams changes with the phase shifts, resulting in changes in fluorescence intensity, which is used as the feedback signal. The phase mask can then be calculated and optimized after a few iterations. Even though these techniques achieve deep tissue imaging, they still face some limitations. First, the maximum improvement in fluorescence intensity for brain imaging is usually less than tenfold [23–25]. Second, the field-of-view after scattering correction is usually limited to a few tens of microns in brain imaging due to the limited range of the memory effect (i.e., isoplanatic patch, mean free path). Third, the total correction speed is relatively slow, including not only the time for taking measurements but also the time for computing the mask, projecting the mask, and reloading updated phase patterns to the SLM for iterative optimization. A detailed comparison between these methods and our technique is summarized in Table 1 (Appendix E).

Active correction of scattering via intensity modulation using a digital micromirror device (DMD) enables high-speed scattering correction. “Intensity modulation” here refers to techniques that directly modulate the intensity of light at the Fourier plane as opposed to “phase modulation” discussed above, which directly modulates the phase of light at the Fourier plane. Most commercially available DMDs have patterning speeds ranging from tens to over 30 kHz, which is significantly faster than the state-of-the-art SLMs. DMDs are also generally more cost-effective than SLMs. Intensity modulation for scattering correction has been demonstrated by focusing a laser beam through a scattering medium and detecting the transmitted laser light with a detector [26–32]. The correction mask can be probed by turning on segments on the intensity modulator one-by-one [26], by projecting Hadamard patterns to calculate the transmission matrix [27,28], or by projecting random patterns combined with transmission matrix measurement [29] or iterative optimization algorithms [30–32]. Compared to phase modulation, traditional intensity modulation faces a challenge: for a given incident laser power on the modulator, scattering correction increases the signal, but intensity modulation reduces the signal because the total intensity on the detector decreases with the turned-off pixels. This is not an issue when focusing a laser beam through scattering media and detecting the transmitted laser light [26–32], as the photon flux from a laser is high.

However, in fluorescence imaging, especially multiphoton imaging of biological samples, scattering correction by intensity modulation is very challenging because: (1) the number of fluorescent photons under two-photon excitation is much lower than that of a laser beam, and the maximum excitation intensity is constrained by the thermal damage threshold of biological samples; (2) fluorescent photons are collected in reflection mode rather than in transmission mode, while most light is forward-scattered by biological tissues, further reducing the number of photons after scattering; (3) two-photon microscopy images embedded objects within scattering tissue, without any free space between the phantom and the focus, which limits the degrees of freedom for correcting scattering [33,34]; (4) imaging a large field-of-view requires scattering correction beyond the memory effect range, which is more difficult than correcting scattering at a single spot.

Recently, high-power femtosecond lasers (ranging from a few watts to over 100 W) for multiphoton microscopy have become commercially available, making two-photon microscopy with intensity modulation feasible [35,36]. It has been used for aberration correction [35,36], but not yet for scattering correction because the challenges mentioned above remain unsolved. In these works [35,36] for aberration correction, the DMD is used as a wavefront modulator by calculating binary holographic masks using Lee holography [37], with Zernike modes employed as the orthogonal basis to probe aberrations. The correction mask is generated using a modified hill-climbing algorithm similar to the classic algorithm in adaptive optics [38], where the final correction mask is a linear combination of weighted Zernike modes, and the weight coefficient of each Zernike mode is iteratively measured. However, binary intensity modulation generated by Lee holography is less efficient and effective than phase modulation, which has been demonstrated for focusing through scattering media [29,39]. Consequently, while this technique corrects aberrations in two-photon microscopy, fluorescence intensity is only improved by 2.5–4 times when imaging ex vivoDrosophila brains at depths of 100–200 µm. Despite the DMD’s fast operational speed of 22.7 kHz, the iterative hill-climbing algorithm is slow, leading to a total processing time of approximately 6 min. Therefore, instead of adapting orthogonal bases (i.e., Zernike modes) and algorithms (i.e., iterative hill-climbing algorithm) designed for phase modulation, innovative approaches tailored specifically for two-photon microscopy with binary intensity modulation are needed, including new orthogonal bases and algorithms. While the key advantage of DMD-based intensity modulation is its patterning speed, no two-photon microscopy system with binary intensity modulation has yet been able to rapidly and effectively correct scattering in deep tissue.

To overcome these challenges, we have developed a new two-photon microscopy system named 2P-FOCUS (Fourier-domain intensity coupling for scattering correction) to provide effective and fast scattering correction over a large field-of-view. The power loss from intensity modulation is compensated for by increasing the laser power input to the intensity modulator (i.e., a DMD) to maintain the same laser power on the sample before and after modulation. Unlike phase modulation, intensity modulation segments the incident collimated beam into multiple narrow beams, generating multiple-beam interference within the scattering tissue. To select the beams that remain in phase after scattering while blocking those that are out of phase, we use random patterns to probe whether the randomly selected beams are in phase or not after scattering. In-phase beams generate a bright main lobe after interference, resulting in brighter fluorescence compared to non-constructive interference. The final correction mask is generated by summing those random patterns followed by binarization. Random patterns are used instead of patterns at specific frequencies (e.g., Hadamard [27,28] or Zernike patterns [35,36]) because: (1) each random pattern consists of a broad range of spatial frequencies, allowing us to multiplex measurements across different frequencies; (2) random patterns can flexibly control the ratio of turned-on pixels to the total number of pixels (“sparsity”) to improve the correction but Hadamard or Zernike patterns cannot; and (3) random patterns generate a central focus with stable intensity and location, whereas Hadamard patterns produce multiple foci that vary in intensity and position with each pattern. This single-shot algorithm ensures constructive interference without the need for time-consuming iterative optimization. To correct scattering beyond the memory effect range, we synchronize the DMD with a galvo-scanner on its conjugate plane to project different correction masks onto corresponding regions during scanning. This strategy significantly increases the system’s degrees of freedom and effectively controls a 4D light field [40].

We demonstrate the impressive capabilities of 2P-FOCUS through a series of proof-of-principle experiments. These include focusing and imaging through approximately 200 µm thick bone and capturing fluorescence-labeled neurons and blood vessels at depths of up to 510 µm deep in ex vivo mouse brain tissue. Remarkably, 2P-FOCUS achieves over a 36-fold increase in fluorescence intensity compared to standard two-photon microscopy when imaging in the brain tissue. 2P-FOCUS also achieves a 230μm×230μm field-of-view in all experiments—beyond the memory effect range and significantly larger than that of current methods. Furthermore, 2P-FOCUS completes the entire correction process in just a few seconds, including the time for measurements, calculation, and projection of correction masks, far surpassing the speed of similar techniques [35,36]. To our knowledge, 2P-FOCUS is the first two-photon microscopy system to use intensity modulation for rapid and effective scattering correction.

2P-FOCUS corrects scattering in three steps: taking measurements, calculating the correction mask, and projecting the correction mask (Figs. 1(A)–1(C)). First, 2P-FOCUS projects random intensity patterns in the Fourier domain while monitoring the fluorescence intensity excited under this modulation (Fig. 1(B), step 1). In the following experiments, random patterns are generated from a uniform discrete distribution on [0,1] using MATLAB, with the random seed shuffled each time a new set of random patterns is created. Different sets of random patterns, consisting of hundreds to thousands of patterns, are generated for varying sparsity levels and super-pixel sizes. Random intensity patterns have been used as an orthogonal bases in computational optics [41]. Random patterns are orthogonal to each other because each of them consists of a large number of segments (hundreds to thousands), thus satisfying the mathematical criteria of an orthogonal matrix, as PTP=I, where P is the reshaped and normalized random pattern matrix and I is the identity matrix. In other words, the dot product of any two distinct random patterns is approximately zero. Each segment on the DMD consists of n×n pixels, termed “super-pixels”. Random patterns are generated using these super-pixels. If N super-pixels with a radius of W are turned on, the electric field after intensity modulation is EA(ξ,η)=∑i=1Nrect(ξ−ξi2W)rect(η−ηi2W),(1)where rect is the rectangular function describing a turned-on super-pixel. The next step is to calculate the electric field E0 on the surface of a scattering medium, assuming the unevenness of the surface within the beam size is negligible compared to the working distance of the objective lens. If the thickness of the scattering medium is Δz0, E0 is calculated by E0(x,y)=P−Δz0F−1EA(ξ,η),(2)where F denotes Fourier transform and P−Δz0 denotes the linear operator for Fresnel propagation backwards by distance Δz0: P−Δz0=F−1exp[−j2πΔz01/λ2−(ξ2+η2)]·F.(3)

That is, the incident field on the surface of the sample is calculated by first computing the electric field at the focus, and then back-propagating it from the focus to the surface over a distance of Δz0. By combining Eqs. (2) and (3), the electric field on the surface of the scattering medium can be simplified to E0(x,y)=F−1exp[−j2πΔz01/λ2−(ξ2+η2)]·EA(ξ,η).(4)

In matrix form, E0=F−1diag(HΔz0*)EA,(5)where H is the forward Fresnel propagation over distance Δz0. The incident electric field E0 propagates through the scattering medium over a distance Δz0, eventually reaching the reference fluorescence object. The excitation field on the fluorescence object is Eexc=TexcE0,(6)where Texc is the transmission matrix of the scattering medium at the excitation wavelength. The intensity of the excitation field Eexc is Iexc=|Eexc|2.(7)

Since each pattern randomly selects multiple beams, they may or may not constructively interfere at a focus after propagating through the unknown scattering media. If these beams interfere constructively, the center lobe of the interference pattern will excite fluorescence through two-photon absorption if the photon density in the center lobe is high enough. Otherwise, the beams do not generate a bright focus, resulting in low or zero fluorescence intensity. Beams that create a bright focus may also be multiply scattered, as long as their phase difference is zero. In two-photon microscopy, fluorescence intensity is quadratically proportional to the intensity of excitation light in two-photon excitation, represented as $$\begin{gathered} If=Iexc2 \\ =|TexcF−1diag(H*)EA|4. \end{gathered}$$(8)

The excited fluorescence propagates backward to exit the scattering medium. Because the wavelength of the fluorescence differs from that of the excitation light, and because fluorescence is incoherent (or partially coherent) light while the excitation light is coherent, we use an intensity transmission matrix Tf for fluorescence. This matrix differs from the transmission matrix Texc used for excitation light. The signal detected by the PMT represents the total fluorescence intensity exiting the scattering medium, and is expressed as $$\begin{gathered} I=∑TfIf \\ =∑Tf|TexcF−1diag(H*)EA|4 \\ =∑Tf|TexcF−1diag(H*)∑i=1NEi|4, \end{gathered}$$(9)where Ei is the matrix form of the binary random patterns with only one super-pixel turned on. Since computing the correction mask by solving two transmission matrices from a nonlinear equation would be extremely challenging, to enable high-speed active correction, our technique probes the correction mask by simultaneously turning on multiple random pixels. During our experiments, several thousand measurements are taken in this process. Since the total fluorescence intensity is proportional to the fourth power of the electric field of these beams, 2P-FOCUS is highly sensitive to the number of pixels that are turned on, resulting in a highly effective and efficient approach for scattering correction.

The next step involves calculating the correction mask (Fig. 1(B), step 2). The random masks that generate bright fluorescence are selected by thresholding the fluorescence intensity detected by the PMT. Since these beams can create bright center lobes through interference, they have little phase difference. Thus, summing these random patterns generates a new pattern. The multiple beams from this new pattern are coherent and in phase, producing a brighter center lobe. This new pattern serves as the grayscale correction mask.

The final step is generating a binary correction mask for DMD projection by thresholding the intensity of the grayscale correction mask (Fig. 1(B), step 3). To ensure a fair comparison, we maintain the same illumination power on the sample before and after correction, staying below the thermal damage threshold of biological tissues. This is a unique aspect of our method compared to conventional intensity modulation. For example, if a blank mask is used before correction, allowing no light blockage (Fig. 1(A)), we then increase the input power to the DMD when applying correction masks to compensate for the power blocked by turned-off pixels (Fig. 1(C)). Power control can also be achieved by using different thresholds to binarize the grayscale correction mask. Different thresholds result in varying illumination power on the sample, given the same input power to the DMD, which will be discussed in detail in the following sections. Since all beams selected by the grayscale correction mask are coherent and in phase, the subset of these beams after binarization remains coherent and in phase. Their interference produces a brighter center lobe compared to the scenario without scattering correction, not only effectively redistributing part of the illumination power from out-of-phase beams to in-phase beams but also nonlinearly increasing the laser power in the center lobe due to coherent interference.

The experimental setup is depicted in Fig. 1(D) and described in Appendix A. We utilize a femtosecond laser with a 1035 nm wavelength for two-photon excitation. Our setup includes a near-infrared coated DMD with 1280×800 pixels, of which only the central 800×800 pixels were used, and a maximum patterning rate of 12.5 kHz. In the following experiments, the DMD is operated at a 1 kHz patterning rate due to limitations of other electronics. Since the DMD is located on the conjugate plane of the back aperture, it induces spatial dispersion in the femtosecond pulses, resulting in an elliptical focus (Fig. 1(E), left). To correct this spatial dispersion and compress the pulses, we implement a reflective diffraction grating on the conjugate plane of the DMD [42,43]. The parameters of the grating, the 4-f relay lenses (L3, L4), and the incident angles are carefully designed and aligned to fully compensate for the spatial dispersion. It is important to note that the effective groove size of the DMD differs from the micromirror pitch size. After dispersion compensation, the point-spread-function (PSF) becomes a tight circular spot (Fig. 1(E), right). A two-axis scanning mirror system, operating at 1 kHz (the maximum speed of the galvo-scanner), is also positioned on the conjugate plane of the DMD. This system scans the focus across the field-of-view and currently represents the speed bottleneck of the system. The scanning mirrors are synchronized with the DMD, allowing the application of different correction patterns to specific subregions as necessary. The modulated light field is relayed to the back aperture via 4-f systems and then focused by the objective lens onto the scattering sample. The sample is placed on a linear translation stage for capturing z-stack images. Two-photon excited fluorescence is collected in reflection mode by the objective lens and detected by a PMT. The 3D PSF of 2P-FOCUS is 1.7μm×1.7μm×6.6μm, measured with 0.71 μm red fluorescent beads without scattering media (Figs. 1(F) and 1(G)). The resolution is slightly lower than theoretical values because the surface of the DMD induces astigmatism, which could potentially be corrected by adding an SLM to the system.

We first demonstrate 2P-FOCUS experimentally by focusing through a chicken bone onto a homogeneous fluorescence slide (Fig. 2). We adhered a piece of chicken bone without any exogenous fluorescent labels onto a microscope slide coated with a thin layer of red fluorescent paint (Fig. 2(E)). The surface of the bone is uneven and the part we focused through is about 200 µm thick. The bone is highly scattering and no focus is visible before correction at an illumination power of 7 mW (Fig. 2(F)). The mean and standard deviation of intensity of Fig. 2(F) is 0.58±1.01. The maximum value of Fig. 2(F) comes from a hot pixel. The random patterns used in the experiment have a sparsity of 0.4 and a super-pixel radius of 8 pixels (see next paragraph for details). During the correction process, we took measurements with 2500 random patterns (Fig. 2(A)) at an illumination power of 76 mW (the power is only used for taking measurements under random patterns) and recorded the corresponding fluorescence intensity (Fig. 2(B)). The random patterns exhibit broad coverage of spatial frequencies, with consistent frequency coverage across all patterns. They maintain the same power in the DC component to ensure high signal-to-noise ratio measurements at the same location, while the power in the high-frequency components varies slightly to probe the spatial components that propagate through the scattering bone (Figs. 6(A) and 6(B), Appendix E). We then selected the random patterns corresponding to the top 10% brightest fluorescence (above the red line, Fig. 2(B), 250 patterns) and summed them to form the grayscale correction mask (Fig. 2(C)). The grayscale value of the correction mask indicates the number of patterns that turn the super-pixel on. For instance, the peak value of 159 in Fig. 2(C), located on the right bottom side in the Fourier domain, shows that most random patterns (159 out of 250) turned this pixel on to create a bright focus. We then adjusted the half-wave-plate (HWP) to maximize the power input to the DMD, and binarized the grayscale correction mask (Fig. 2(D)) using a threshold that ensures the illumination power on the sample remains at 7 mW after modulation by the binary correction mask, which is the same as before correction. In this case, the correction mask only turns on pixels that can generate highly constructive interference (see next paragraph for details). Given that the turned-on pixels are mostly on the edge of the back aperture in this case, our interpretation is that the bone is probably uneven in the field-of-view and the correction mask tends to guide light through the relatively thinner (less scattering) part. The correction of a relatively flat bone is shown in Fig. 3. This correction mask leads to a bright fluorescent focus with a peak intensity of 629 (Fig. 2(G)), measured with a camera in the reflection mode (replacing the PMT in Fig. 1(D)) for visualization purposes only. The distorted shape of the fluorescence spot is due to the scattering of fluorescence as it propagates back through the bone, as the correction applies only to the excitation beam and not the emitted fluorescence. This is not an issue in two-photon microscopy, as scattered emission fluorescence can still be collected by the PMT as long as it falls within the numerical aperture (NA) of the objective lens. The intensity profile across the x-axis (Fig. 2(H)) also demonstrates the significant improvement in fluorescence intensity before and after correction.

We next explore how various factors influence the correction result, including the thresholds for binarization, the sparsity of random patterns, the radius of super-pixels, and the number of measurements. The binarization in Fig. 2(C) and Fig. 2(D) depends on the threshold applied to the grayscale correction mask, which controls the illumination power applied to the sample given the same input power to the DMD. The experiment depicted in Figs. 2(A)–2(H) describes a scenario where the laser power on the sample is 7 mW, corresponding to the first data point in Fig. 2(I). The fluorescence intensity value shown in Fig. 2(I) is measured by the PMT, so it differs from the absolute values in Figs. 2(F)–2(H). The “before” group in Fig. 2(I) is measured under random patterns, and the power for the “before” group is adjusted by tuning the HWP. Random patterns instead of a blank screen are used here to check whether the random pattern itself causes a non-quadratic relationship between the illumination power and the fluorescence intensity. As the illumination power increases (to 15 mW, 34 mW, 52 mW, and 76 mW), the fluorescence intensity before correction increases quadratically with the illumination power (blue dashed line in Fig. 2(I), and Fig. 6(C)), consistent with the principles of two-photon excitation. Importantly, after applying the correction mask, the intensity of two-photon excited fluorescence does not follow the quadratic rule as the excitation power increases by turning on more pixels on the DMD (blue solid line in Fig. 2(I), and the corresponding correction masks are shown in Figs. 6(D) and 6(E)). This means the significance of each pixel on the DMD is different. The brightest pixels in Fig. 2(C) generate “most in-phase beams”, while other bright pixels generate “partially in-phase beams”. Turning on the DMD pixels corresponding to the partially in-phase beams increases the illumination power on the scattering sample, but only part of the laser power is delivered to the fluorescence object and constructively interferes with other beams. In addition to partially in-phase beams, several other factors can also increase the power on the sample without contributing to two-photon fluorescence. One such factor is the activation of pixels corresponding to low-transmission modes, as each pixel on the DMD corresponds to different spatial modes with varying transmission through the tissue (Fig. 2(C)). These beams may also be diffracted away from the focus of constructive interference. The most significant improvement in fluorescence intensity is achieved when the laser power on the sample is relatively low. In other words, the strategy to achieve a large improvement after correction is making the correction mask turn on only the pixels generating highly in-phase beams after scattering (red line, Fig. 2(I)). This result clearly proves the mechanism of 2P-FOCUS is selecting the multiple beams that are still highly in phase after scattering while blocking those out-of-phase ones.

The second factor is the sparsity of random patterns, which indicates the percentage of pixels that are turned on out of the total number. For random patterns with a super-pixel radius of 8 pixels, there are 50×50=2500 super-pixels in each mask, and a total of 2500 masks are used, forming a 2500×2500 matrix P. The sparsity of random patterns ranges from 0.05 to 0.9 (see examples in Fig. 7(A), Appendix E). To demonstrate that the random patterns are orthogonal to each other, we calculated I=PTP using normalized random patterns with a sparsity of 0.4. The histogram of all elements in I shows that all diagonal values equal one and non-diagonal values are 0±0.02, indicating that I is approximately an identity matrix (Fig. 7(B)). Consequently, the random patterns are orthogonal to each other. The sparsity of random patterns is also related to illumination power. With the same input laser power to the DMD, less sparse patterns block less input power and emit more power to the sample. This increases the SNR of the measurements, consequently improving the accuracy of the correction masks. In the case shown in Fig. 2(J) (the values are listed in Table 2, Appendix E), we first measure the fluorescence intensity using random patterns at 7 mW, as in the pre-correction case, and then take measurements under different sparsities (0.05, 0.1, 0.2, 0.3, and 0.4), while keeping the input power to the DMD constant. These measurements generate five different correction masks, each with an output power of 7 mW (see Fig. 7(C)). We also took measurements under sparsities ranging from 0.2 to 0.5 through another region of the bone and found that the fluorescence signal decreases at the sparsity of 0.5 (see Figs. 7(D) and 7(E) and Table 3, Appendix E). The correction mask derived from the random patterns with 0.4 sparsity shows the most significant improvement in fluorescence intensity in both cases, as depicted in Figs. 2(J) and 7(D). In practice, the optimal sparsity of random patterns varies depending on the scattering properties of the sample, the brightness of the fluorescence signal, and the maximum power of the laser.

The third factor influencing the correction result is the size of super-pixels. The center square region of the DMD consists of 800×800=640,000 pixels. Well-sampling all these pixels requires 640,000 random patterns, which is too large for the on-chip RAM of the DMD to handle and also requires more computing and experimental time. Therefore, we binned these pixels into super-pixels to accelerate the process. For super-pixels with a radius of W, the number of random patterns to well-sample the frequency domain is 8002/(2W+1)2. We compared the results generated by patterns with W=8, 12, and 21 (see Fig. 2(K)). The fluorescence intensity before correction was measured under random pattern modulation with a sparsity of 0.05. The illumination power on the sample was approximately 12 mW, similar before and after correction for all three pattern types. Before correction, fluorescence intensity increases as the super-pixel size increases due to the effect of diffraction. Focusing a beam modulated by random patterns on the back aperture generates a main lobe of focus with speckles due to the diffraction of the random patterns. Random patterns with smaller super-pixel sizes diffract more than those with larger super-pixel sizes, resulting in lower peak power at the main lobe and higher power at the speckles given the same total input power. After correction, correction masks with smaller super-pixels correct scattering more effectively because they identify in-phase beams with greater precision.

We also explored how the number of measurements influences scattering correction. For random patterns with a super-pixel radius of 12 pixels, the minimum number of measurements required for well-sampling is 1024. We generated correction masks at 0.01 sparsity with various numbers of raw measurements, ranging from 200 to 2000, and compared the results of the scattering correction. As shown in Fig. 2(L), when the number of measurements is too small (i.e., 200 measurements in this case) to reconstruct an accurate correction mask, this failed mask generates a main lobe of focus and speckles on the image plane, similar to a random pattern. This leads to the fluorescence excited by the focus being even lower than that without correction, which projects a blank mask on the back aperture like standard two-photon microscopy. The blank mask is used to measure the baseline, as it is commonly used in standard two-photon microscopy. The fluorescence intensity improves with an increasing number of measurements, but this improvement becomes marginal when the measurements become redundant. Considering the computing and experimental time, the results suggest that the optimal number of measurements corresponds to the minimum needed to adequately sample the Fourier domain.

We next explore the optimal strategy for imaging a large field-of-view through highly scattering media using 2P-FOCUS. 2P-FOCUS can correct scattering during the imaging process with point scanning; it scans the focus across the field-of-view pixel by pixel while applying correction masks. The most straightforward strategy is to apply the same correction mask for all scanning positions, referred to as “global correction” (Fig. 3(A)). Global correction is valid if the scanning area is within the memory effect range [44]. However, for severe scattering media, the scattering properties are highly spatially varying, resulting in a limited memory effect range. The range is typically a few tens of microns for highly scattering biological tissue [13], depending on the thickness of the scattering layer and the wavelength of the incident light. In this case, the field-of-view for imaging through scattering tissue is limited by the memory effect range.

To image beyond the memory effect range, we developed a new strategy called “subregion correction” (Fig. 3(B)), which synchronizes the scanner and the DMD to seamlessly project different masks to different regions during the scanning process. To ensure the speed of correction, we developed a control code that can automatically perform subregion correction. First, it captures the image before correction, splits the whole image into multiple subregions (e.g., 3×3), and automatically selects the pixel with the maximum intensity as the reference for each subregion, which can be fluorescence-labeled objects or background fluorescence. If the fluorescence intensity of a reference is too low, indicating it may be noise, the subregion is merged into the adjacent one. Next, the code automatically directs the scanner to each reference and generates the corresponding correction mask. The same set of random patterns is used in every subregion and preloaded into the DMD’s RAM to minimize data communication time. After generating all subregion correction masks, the code synchronizes the scanner and the DMD to apply each subregion correction mask to the corresponding region while the sample remains steady. The scanner operates at its maximum speed of 1 kHz during the process. This seamless and automated control significantly accelerates the subregion correction process, making scattering correction beyond the memory effect range feasible in practice.

We experimentally compared global correction and subregion correction side-by-side by imaging red fluorescent beads through a piece of chicken bone (Fig. 3(C)). The size of the fluorescent beads is 0.71 µm, suspending in PDMS, with some clusters. The bone we imaged through is about 200 µm thick. The field-of-view of imaging is 230μm×230μm, and the post-objective power is 25 mW for all image conditions. We first imaged the sample without any correction by projecting a blank mask at all scanning locations. The image (Fig. 3(E)) shows a low signal-to-noise ratio, revealing two bead clusters and several single beads with severe aberrations. Because the bone on top is highly heterogeneous, the aberration and scattering caused by the bone are also spatially varying. We next performed a global correction using the bead cluster indicated by the yellow arrow in Fig. 3(E). The correction used 10,000 random patterns with a super-pixel size of 4 pixels and a sparsity of 0.4. By applying the same correction mask to all scanning locations, the fluorescence intensity after global correction improved by 2.9 times (Fig. 3(F)). However, only the fluorescent bead cluster in the area where the correction mask was measured becomes brighter (indicated by the yellow arrow in Figs. 3(E) and 3(F)), while the fluorescent bead cluster that was 96 µm away (indicated by the red arrow in Figs. 3(E) and 3(F)) does not.

Next, we conduct subregion correction on the same field-of-view. The image taken under global correction is divided into 3×3 subregions (Fig. 3(D)), and the peak intensity of each subregion is automatically selected as the reference for that subregion (pointed out by yellow arrows in Fig. 3(D)). Both fluorescent beads and background fluorescence can serve as references to generate correction masks, with the selection being automatically made by the code based on intensity. Brighter objects provide higher signal-to-noise ratio measurements, resulting in more accurate correction masks. The same set of random patterns used in global correction is applied to generate the correction mask for each subregion. Under subregion correction, the image (Fig. 3(G)) shows improved intensity and enhanced contrast across the entire 230μm×230μm field-of-view. Objects that were previously invisible under global correction become visible (indicated by the red arrow in Figs. 3(E)–3(G)). The peak intensity improved by 4.8 times compared to the pre-correction case. The intensity profile across the two bead clusters quantitatively compares the results of global and subregion correction (Fig. 3(H)). The line plot indicates that the memory effect range of this bone sample is less than 96 µm, while the thickness of the bone is about 200 µm. Therefore our method corrects multiple scattering. Additionally, to examine potential photobleaching due to the correction process, we imaged the same region with a blank mask before and after the correction process (Fig. 8, Appendix E). The result shows the peak intensity slightly decreased from nine to eight, indicating only minor photobleaching. Both global and subregion corrections effectively improve the fluorescence intensity of objects within the memory effect range (left peak in Fig. 3(H)), while only subregion correction improves intensity beyond the memory effect range (right peak in Fig. 3(H)).

We next check the image resolution using a single bead before correction, after global correction, and after subregion correction (Figs. 3(I)–3(K)). Before scattering correction, the image of a single bead is elongated due to strong aberrations from the bone, with an FWHM of 2.48 µm along the short axis and an average peak intensity of 4.9 (Fig. 3(I)). After global scattering correction, where the reference is outside of the memory effect range, the bead still looks elongated with an FWHM of 1.91 µm and an average peak intensity of 8.23 (Fig. 3(J)). The slight increase in fluorescence intensity is probably because the scattering potential in this region is coincidentally similar to the region where the correction is conducted. After subregion correction, the bead is still elongated with an FWHM of 1.87 µm, but the peak intensity improves to 31.55 (Fig. 3(K)). Unlike phase modulation or intensity modulation for aberration correction [35,36], our method selects in-phase beams to maximize fluorescent intensity passing through scattering media without correcting aberrations. Our technique could be further enhanced when combined with aberration correction methods, which is out of the scope of this paper. The results quantitatively demonstrate the intensity modulation does not reduce spatial resolution compared to the case before correction when all pixels on the DMD are turned on, because the blocked beams cannot pass through the scattering sample and constructively interfere even if the corresponding DMD pixels are turned on.

To further demonstrate our approach for correcting spatially varying scattering, we use subregion masks generated with different references (Fig. 3(D)) as global correction masks to scan across the entire field-of-view. The “on” pixels in the correction patterns (Fig. 3(L)) are predominantly located near the center, likely due to the Gaussian beam profile and global aberration. Assuming the surface of the sample is relatively flat, on-axis ballistic photons propagate a shorter distance than off-axis ballistic photons, resulting in more on-axis ballistic photons reaching the focus. On the other hand, the correction masks function differently in correcting scattering at various locations (Fig. 3(M)) because of their high-frequency features. For example, mask #1 and mask #4, calculated using the left bead cluster (pointed out by the yellow arrow in Figs. 3(E)–3(G)), brighten the left cluster but not the right, while mask #7 and mask #8, calculated with the right bead cluster (pointed out by the red arrow in Figs. 3(E)–3(G)), do the opposite. Masks calculated from relatively dim references, like #2, #3, and #9, show less effective correction. The results indicate that 2P-FOCUS improves fluorescence intensity by scattering correction, since the mask tailored to one specific scattering region cannot correct scattering in another region beyond the memory effect range. Subregion correction is necessary to address scattering outside the memory effect range across a large field-of-view.

To further explore the contribution of high-frequency and low-frequency components in the correction masks, we applied various smoothing filters to the correction masks to partially remove the high-frequency components (Fig. 8(A), plots 1–5). Additionally, we applied a low-NA mask as a limiting case that excludes all high-frequency components (Fig. 8(A), plot 6), which has been used to enhance imaging depth in two-photon microscopy [45,46]. The illumination power on the sample was kept constant under all masks. The corresponding fluorescence intensity under each mask (Figs. 8(B) and 8(C)) demonstrates that maximum fluorescence intensity is achieved after slight smoothing (plot 4). In contrast, both the absence of smoothing and excessive smoothing (plots 5 and 6) result in a drop in fluorescence intensity. Without smoothing, the high-frequency components diffract light severely, reducing the power in the central lobe. On the other hand, excessive smoothing or the low-NA mask fails to effectively correct highly spatially varying scattering, also reducing the power of the focus. Plot 3 and plot 4 exhibit similar intensities, indicating that the positions of some super-pixels are less significant than others. This observation aligns with our findings in Section 3.A, where we noted that beams from certain super-pixels can constructively interfere while others cannot. These results demonstrate that 2P-FOCUS primarily improves fluorescence intensity by correcting highly spatially varying scattering, with aberration correction playing a secondary role, and fundamentally differs from approaches that enhance imaging depth using a low NA.

We next demonstrate the imaging capabilities of 2P-FOCUS by imaging fluorescence-labeled neurons in the mouse brain. In our first imaging experiment, we imaged a 230μm×230μm×450μm volume in the primary visual cortex of an intact transgenic mouse brain with cell-fill expression of tdTomato in all parvalbumin (PV) interneurons (see Fig. 4(A)). The mouse underwent transcardial perfusion with 4% paraformaldehyde (PFA) and was then preserved in 4% PFA until imaging (see Appendix D for detailed sample preparation). During the imaging session, the whole brain was maintained in phosphate-buffered saline and placed beneath the objective lens. The brain was imaged from the top downward. The top plane (z=0) was not on the brain’s surface but below it, as dendrites and axons are observed on the plane (Fig. 4(B)). When imaging the top 250 µm, we did not apply any correction but gradually increased the illumination power on the sample from 2 mW to 14 mW. The first correction was applied when imaging the plane at a 250 µm depth (see Figs. 4(D)–4(F)). Before correction, the image contrast is low, but we can still identify some bright neurons (Fig. 4(E)). We selected one of the bright neurons near the center of the field-of-view, indicated by the yellow arrow in Fig. 4(E), as the reference object. It is important to note that the reference neuron is not a guide star, as it does not need to be a point source; any fluorescent object can be used as a reference. We applied an offset voltage to the scanning mirrors to focus on the reference neuron, and then performed a scattering correction to maximize the fluorescence intensity at this neuron. We used 2500 random masks with a super-pixel radius of 8 pixels and a sparsity of 0.4 for measurements, and generated a correction mask as shown in the bottom of Fig. 4(D). We applied this correction mask to all scanning positions as a global correction, and acquired the image of the same region after correction (Fig. 4(F)). The maximum fluorescence intensity is improved by about 3-fold (from 27 to 80). The intensity of the reference neuron is improved by about 1.8-fold (see Figs. 9(A)–9(C), Appendix E), which is less than the maximum improvement and may indicate minor photobleaching during the correction process. The image contrast is also improved after correction, as the three neurons in the zoomed-in view (Figs. 4(E) and 4(F)) were not distinguishable before correction but became distinguishable after correction. The intensity profiles of the three neurons are plotted in Fig. 4(K).

As the imaging depth increases, the global correction mask no longer provides high-intensity and high-contrast images when the imaging region is larger than the range of the memory effect. We applied subregion correction to image deep in the mouse brain. We first conduct subregion correction at 290 μm where the peak intensity drops to the half of the peak intensity of the image at 250 μm depth after global correction (see Figs. 9(D)–9(F)), indicating the image region at 290 μm depth is out of the memory effect range from the image region at 250 μm. The field-of-view is divided into three subregions, and three correction masks are generated using reference fluorescence at 290 μm (see Fig. 9(G)). Notice that one of the reference fluorescence is the background fluorescence rather than a neuron. The three correction masks are applied when imaging the next few z planes (z=290–380μm). The peak intensity of the image at 380 μm using the correction masks generated at z=290μm is 66 (see Fig. 9(J)), indicating this imaging area was out of the memory effect range from the imaging area at 290 μm depth. Thus, we conduct another round of subregion correction at 380 μm depth and generated four correction masks (Figs. 9(J)–9(L)). These masks are then applied for imaging z planes at 380–450 μm depth. Figures 4(H) and 4(I) show the image at 400 μm depth before and after subregion correction. The subregion correction successfully improved the fluorescence intensity by 9.3 times (see Fig. 4(I)) compared to the image without correction (Fig. 4(H)). The neuron barely visible before correction becomes visible, indicating an improvement in the SNR (see zoomed-in views in Figs. 4(H) and 4(I)). The intensity profile of the representative neuron pointed out by the arrows is shown in Fig. 4(L). The improvement is less significant than that achieved when focusing through bone, because the correction mask used for focusing through bone is more selective than the one used for imaging neurons. The fluorescence slide used in the focusing experiment is much brighter than the fluorescence-labeled neurons in the mouse brain. Given the same input laser power, with brighter fluorophores, the correction mask can selectively turn on only the pixels corresponding to highly in-phase beams and still achieve a sufficient SNR. In contrast, with less bright fluorophores, the correction mask has to turn on more pixels, including those corresponding to partially in-phase beams. Therefore, the intensity improvement in brain imaging is less than that in focusing through bone, which is consistent with our result in Fig. 2(I). The bottleneck is the power density of the input laser to the DMD. Subregion correction enables us to image a 230μm×230μm region deep within the scattering brain tissue, achieving better results than global correction.

To further improve the correction and demonstrate the broad applicability of 2P-FOCUS, we used finer random masks with a super-pixel radius of 4 pixels at 0.4 sparsity and imaged another mouse brain ex vivo, in which cerebral blood vessels are labeled via intravascular FITC-dextran injection (Fig. 5(A)). Similar to the setup in Fig. 4, the whole mouse brain was immersed in phosphate-buffered saline and imaged from the top downward. The top plane was not on the brain’s surface but below it, as indicated by the visible cross-sections of cerebral blood vessels. The top 270 μm volume was imaged without correction, meaning all pixels of the DMD are turned on. Since blood vessels span three-dimensionally, we choose to show the maximum intensity projection rather than a single plane. The maximum intensity projection of the top 100 µm thick volume is shown in Fig. 5(B), where the cerebral blood vessels are clearly visible. As the imaging depth increases, the laser power on the sample is gradually increased from 3 mW (for z=0–160μm) to 19 mW (for z=250–510μm). The first scattering correction is conducted at a depth of 270 μm in three subregions, using 2500 random patterns with a super-pixel radius of 8 pixels in each subregion. At a depth of 300 μm, we conducted another subregion correction with four subregions, using 8100 random patterns with a super-pixel radius of 4 pixels in each subregion. Compared to the image before correction (Fig. 5(D)), the image after correction (Fig. 5(E)) reveals fine structures of cerebral capillaries that are invisible before correction. The peak intensity is improved by a factor of 36.1 after correction, which is more than that of the state-of-the-art techniques for correcting scattering in the mouse brain tissue (Table 1). The four correction masks and subregions are shown in Figs. 5(I) and 5(J). The same masks were used when imaging at a depth of 340 μm (Figs. 5(F) and 5(G)). Compared to the uncorrected image with the same illumination power (Fig. 5(F)), where all pixels on the DMD were turned on, these capillaries were barely visible. The fluorescence intensity of a cross-section of representative capillaries improves by a factor of 30.6 after correction (Fig. 5(H)). Another example for scattering correction at 390 μm depth with six subregions is shown in Figs. 10(A)–10(C) (Appendix E). At the bottom of the volume (410–510 μm depth), blood vessels remain visible after scattering correction (Fig. 5(C)). Image resolution shows no obvious change after applying the correction mask, as fine structures like capillaries with a diameter of around 2 μm are resolvable (Figs. 10(D)–10(H)). Compared to the results in Fig. 4, the outcomes in Fig. 5 show better scattering correction, attributable to the use of patterns with smaller super-pixels that generate more precise correction masks. Additionally, the fluorescence labeling of blood vessels is brighter than that of the neurons, despite 1035 nm not being the optimal two-photon excitation wavelength for FITC-dextran. Furthermore, in the experiment, the mouse brain is not fixed by 4% PFA, which results in less autofluorescence compared to the fixed mouse brain shown in Fig. 4. Therefore, the SNR of measurements under random pattern modulation is higher, resulting in a more accurate correction mask.

We have demonstrated that 2P-FOCUS is a novel and fast method for correcting scattering in two-photon microscopy. The key advantages of 2P-FOCUS are its speed and effectiveness, achieved through the innovative joint design of hardware and software. In terms of hardware, 2P-FOCUS uses a DMD, which has a much faster patterning rate than an SLM, to select in-phase beams, enabling high-speed patterning at a lower cost. In terms of software, 2P-FOCUS utilizes random patterns as an orthogonal basis and a single-shot algorithm to complete the entire correction process within a few seconds depending on the number of measurements. For example, collecting 2500 measurements at a 1 kHz sampling rate takes 2.5 s, calculating the correction mask takes 0.5 s, and projecting the correction mask takes 1 ms—a total of 3 s to complete the entire correction process, with no iterations needed. In each subregion correction task, the random patterns used for all subregions are identical and preloaded into the DMD’s RAM, which significantly reduces the time for data transfer. Different sets of random patterns are employed across various experiments to minimize the likelihood of accidentally favorable patterns and to demonstrate the robustness and consistency of 2P-FOCUS. 2P-FOCUS effectively improves signal intensity by leveraging the nonlinearity of multiple-beam interference and two-photon excitation. Unlike phase modulation, intensity modulation selects multiple in-phase beams that constructively interfere after passing through scattering media. The nonlinearity of multiple-beam interference greatly enhances the intensity in the central lobe compared to the linear superposition of multiple incoherent beams. Two-photon excitation further amplifies this nonlinearity compared to one-photon excitation. Overall, the fluorescence intensity is proportional to the fourth power of the multiple-beam field in the Fourier domain, making the fluorescence intensity extremely sensitive to the intensity modulation. Although modulating the intensity of light at the Fourier plane can induce a phase change at the image plane (as in Lee holography), 2P-FOCUS calculates the intensity mask directly from the fluorescence intensity of multiple-beam interference rather than relying on phase correction using Zernike modes. This approach classifies our method as an intensity modulation technique rather than a phase modulation technique. 2P-FOCUS also enables subregion correction and highlights the importance of subregion correction for addressing scattering beyond the memory effect range. We demonstrated 2P-FOCUS by correcting scattering in mouse brain tissue ex vivo over a 230μm×230μm×510μm volume, successfully improving the fluorescent intensity by 36.1-fold after correction.

To apply 2P-FOCUS for in vivo imaging in the future, the current 2P-FOCUS needs to be improved in two key aspects. First, to extend the imaging depth further, 2P-FOCUS needs to acquire high signal-to-noise ratio measurements even deep inside the sample. This proves challenging for the current version of 2P-FOCUS, as measurements are taken under completely random modulation in the Fourier domain. Sufficient signal-to-noise ratio of measurements is necessary for calculating an accurate correction mask, as shown in Fig. 2(J). This limitation could potentially be overcome by incorporating the correction mask from the plane above to generate the correction mask for the current plane. Second, the current process of generating a correction mask in 2P-FOCUS is not fast enough to correct millisecond-scale dynamic scattering due to blood flow [6]. The current correction speed may be slow for imaging through blood flow, but it should be sufficient for imaging areas without arteries, as previous works [23] based on wavefront shaping show that the correction mask can last for tens of minutes in vivo. The maximum patterning rate of the DMD is 12.5 kHz, which means the correction process could take less than 1 s if the DMD could be operated at full speed with upgraded electronics. Therefore, the correction speed of 2P-FOCUS could be further improved with a higher bandwidth data acquisition card, a higher repetition rate laser, and software designed for auto-correction with closed-loop control.

In summary, 2P-FOCUS is the first two-photon microscopy system to use intensity modulation for rapid and effective scattering correction. It also uses random patterns as the orthogonal basis and employs a single-shot algorithm to complete the correction within a few seconds. This method can significantly boost fluorescence intensity by several tens of times at the same illumination power on the sample. This capability makes it a powerful tool for deep tissue imaging, facilitating research in neuroscience, immunology, and cancer.

Acknowledgment. We thank Hillel Adesnik and Janine Beyer at the University of California, Berkeley, for preparing the transgenic mouse brain in Fig. 4 and for providing the PMT used in the microscopy setup. We thank Jiandi Wan and Brianna Urbina at UC Davis for preparing the mouse brain in Fig. 5. We thank Weijian Yang and Zixiao Zhang at UC Davis for helping us measure the pulse width of the laser light before and after dispersion compensation. We thank Peter T. C. So at the Massachusetts Institute of Technology and Laura Marcu at the University of California, Davis, for proofreading the manuscript. Research reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. This work was also supported by Dr. Yi Xue’s startup funds from the Department of Biomedical Engineering at the University of California, Davis.