Dual deconvolution in multiphoton structured illumination microscopy for deep-tissue super-resolution imaging

Domain: chemistry
Published: 2026-03-04
Editorial outcome: major_revision

Abstract
Imaging in thick biological tissues is often degraded by sample-induced aberrations, which reduce resolution and contrast, particularly in super-resolution techniques. While hardware-based adaptive optics (AO) using wavefront shaping can correct these aberrations, their complexity and cost hinder widespread adoption. Here, we present a computational AO framework for multiphoton structured illumination microscopy, enabling deep-tissue super-resolution imaging with minimal hardware modifications. By replacing the photodetector with a camera from the conventional laser-scanning multiphoton microscope, we capture a sequence of scanned images. Using virtual structured illumination, we develop a dual deconvolution algorithm that independently corrects excitation and emission aberrations, recovering an aberration-free object spectrum with an extended spatial frequency bandwidth. We experimentally validate this framework through two-photon super-resolution imaging, achieving a lateral resolution of 130 nm—one-fourth of the emission wavelength—at a depth of 180 μm in thick mouse brain tissue, where conventional deconvolution fails to maintain super-resolution capability. This approach provides a cost-effective and accessible alternative to hardware-based AO, expanding the potential for high-resolution deep-tissue imaging in biological research.

Introduction
Fluorescence imaging is an essential tool in biological research, offering molecular specificity and high contrast. Over the past few decades, significant efforts have been made to overcome the diffraction limit of resolution, leading to the development of various super-resolution (SR) imaging techniques. Stimulated emission depletion (STED) microscopy achieves resolution enhancement by exploiting nonlinear point spread function (PSF) engineering1, while single-molecule localization microscopy (SMLM) reconstructs high-resolution images by precisely localizing individual fluorescent molecules2,3. Super-resolution optical fluctuation imaging (SOFI) improves resolution by analyzing higher-order temporal correlations in fluorescence signals4. Structured illumination microscopy (SIM) doubles the diffraction-limited resolution using patterned excitation, offering ~150 nm lateral and ~300 nm axial resolution5. Although SIM has lower resolution than STED or SMLM, it provides advantages in imaging speed and compatibility with conventional fluorophores.

A key challenge in SR imaging is extending the imaging depth to visualize fine structures within biological tissue6. Conventional SIM is generally limited to shallow depths of 10–20 μm in mouse brain tissue. To improve penetration depth, point-scanning SIM techniques such as multifocal SIM (MSIM)7 and instant SIM (iSIM)8 have been developed, enabling imaging depths of ~50 μm in live zebrafish embryos9. By integrating two-photon excitation with line-scanning SIM10 or point-scanning SIM, imaging depths exceeding 100 μm have been achieved in the eye of a live zebrafish embryo8. However, tissue-induced aberrations distort the PSF, degrading both resolution and imaging depth. This issue is particularly critical for SR techniques, which are highly sensitive to PSF distortions11,12.

To resolve tissue-induced aberrations and further extend SR imaging depth, hardware-based adaptive optics (AO) has been incorporated into SIM13–16. AO uses wavefront control devices, such as deformable mirrors or spatial light modulators, to correct aberrations in the excitation path16,17. These aberrations can be measured using direct wavefront sensing (e.g., Shack-Hartmann sensors)17,18 or sensorless approaches13 that iteratively optimize image metrics. Hardware-based AO-SIM has demonstrated imaging depths of ~25 μm in mouse brain tissue and ~100 μm in transparent zebrafish larvae13. When combined with two-photon excitation, AO extends imaging depths up to ~250 μm14. However, because excitation and emission wavelengths differ, they experience distinct aberrations, requiring separate corrections. Recent advancements in two-photon AO-SIM have enabled dual-path (both excitation and emission path) aberration correction, achieving imaging depths up to 500 μm in mouse brain tissue19. Despite these improvements, hardware-based AO systems introduce significant complexity, require additional calibration steps, and remain inaccessible to many biological researchers due to cost and technical constraints.

Computational AO provides an alternative by correcting aberrations directly from conventionally acquired fluorescence images, without requiring additional optical hardware. However, fluorescence imaging lacks coherent wavefront information, making it challenging to directly retrieve aberrations as in coherent imaging techniques20–22. Image reconstruction methods, such as deconvolution23–28 or machine-learning-based AO method29,30 are commonly used to estimate the most likely PSF and sample image from a single fluorescence image. Recently, computational aberration compensation on the detection path of conventional wide-field fluorescence imaging has been demonstrated by constructing a fluorescence reflection matrix31, the incoherent analog of the coherent reflection matrix32. Nevertheless, these methods do not independently account for the excitation and emission PSFs. Instead, they either consider only the emission PSF or adopt a single effective PSF17, which is defined as the product of excitation and emission PSFs. This results in a narrower spectral frequency bandwidth than that of individual PSFs, leading to the loss of high-frequency information. As a result, conventional deconvolution methods perform well under mild aberrations but struggle in conditions with strong aberrations.

Here, we introduce a computational AO framework for multiphoton SR imaging that overcomes these limitations. Our approach builds upon image scanning microscopy (ISM)33, replacing the conventional integral detector with a camera in a laser-scanning microscope. On the basis of multiphoton virtual structured illumination31, we develop a dual deconvolution algorithm that independently corrects the excitation and emission PSFs. By separately accounting for these two sources of aberration, we reconstruct the aberration-free object spectrum with an extended spatial frequency bandwidth. Unlike conventional deconvolution methods that use vector decomposition to find a single effective PSF, our dual deconvolution technique applies matrix decomposition to correct excitation and emission PSFs independently. This significantly enhances image reconstruction, particularly for high-frequency components, and enables SR imaging even under severe aberrations. We validate our framework in two-photon fluorescence microscopy, demonstrating the visualization of dendritic spines in mouse brain tissue with a lateral resolution of 130 nm, corresponding to a quarter of the emission wavelength, at an imaging depth of 180 μm. By providing a cost-effective and hardware-free alternative to AO, our approach has the potential to broaden access to deep-tissue SR imaging for biological research.

Results
Multiphoton structured illumination imaging in the aberrating medium
The proposed method begins with an imaging configuration where the sample is illuminated with a focused beam

Methods

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This article by Choi and Yoon presents a unique approach to achieve super-resolution imaging in multiphoton scanning microscopy in aberrating media. The authors use a modified version of a technique known as image scanning microscopy (ISM) that captures images of the emitted point spread function in a laser scanning microscope setup. They demonstrate their novel approach of leveraging the diversity of illumination intensity angles and phase shifts present in the imaged PSF that naturally emerge during point scanning to perform an alternate version of structured illumination. This enables dual-deconvolution computational adaptive optics for the excitation and emission light. Their calculated theoretical resolution is 140nm, which is very impressive for distorting media. In practice, they find closer to 175nm through a scattering phantom and 240nm without scattering media. The manuscript is thorough, and presents a compelling case for the novelty and impact of the proposed concept. However, some additional experiments are needed to conclusively demonstrate the technique and to compare with similar current best practices. One notable lack is in the comparison between the described approach and standard ISM. Further discussion of this would help to strengthen the manuscript. Although the manuscript requires revision, the approach is cutting edge and constitutes a significant advance in the use of computational and hardware approaches to create super-resolution images in living animals. It seems to be of sufficient impact for publication in Nat Commun and is likely to be of interest to the journal’s readership.\nAdditional comments:\n-Given that this paper presents an enhanced version of image scanning microscopy, a direct comparison in the theoretical and attainable resolution versus typical ISM pixel reassignment reconstructions that have been implemented in the past would be relevant to include in the text and figures.\n-Lines 122-125— While the details of the dual-deconvolution in the methods seems to make sense to this reviewer, the description in the results is harder to follow. It is seems reasonable that the idea of tuning the phase of each spectrum to enhance the intensity of the synthesized spectrum will help with the excitation aberrations. Figure 1 does a decent job of alluding to this possibility. However, it is somewhat less clear in the text and figure 1 how the modification of the spectral map (as described in lines 137-138) is performed for the emission.\n-Lines 160 and following, 403 and following, and Figure 2— The authors state that “the excitation and emission pupil maps were computationally prepared by the random superposition of Zernike modes numbering up to 30.” How many random combinations were tested? What type/depth of tissue aberrations is this hoping to mimic? Were the random aberrations different between emission and excitation, and how realistic is this combination?\n-Are signal loss, noise and tissue scattering (or even aberrations beyond Zernike mode 30) likely to play a role in the limitations of the approach? This is mentioned briefly in the discussion, but additional information could be helpful.\n-Figure 3 and lines 208-219— The bead separation measurement is insightful, but what was the experimentally measured size of the sub-diffraction beads? How close does this come to the theoretical predictions?\n-It would be interesting to see the calculated excitation and emission Zernike modes for some of the samples displayed. How stable are the optimized Zernike modes for a given set of aberrated images?\n-Figure S3— Since the idea of isoplanatic patch is mentioned, it could support the concept to calculate the isoplanatic patch size for the image shown. This could be accomplished by applying the optimized Zernike modes to neighboring locations.\n-line 312— A little bit more discussion about other approaches to speed imaging in ISM and alternate detectors (e.g. SPAD cameras) prior to the concluding paragraph could help to strengthen the case for future impact of the work.", "reviewer_2": "This manuscript presents a solid advancement in multiphoton super-resolution imaging through the development of a computational adaptive optics (AO) framework based on virtual structured illumination microscopy (SIM). A central innovation is the dual deconvolution algorithm, which iteratively corrects excitation and emission aberrations using Wiener-like updates to solve for three unknowns (i.e., object, excitation OTF, and emission OTF), thus decomposing a complex optimization problem into well-conditioned, tractable subproblems. While the spectral synthesis concept underlying the two-photon SIM implementation is similar to the recently established k-space interpretation of image scanning microscopy in the single-photon regime, the proposed dual deconvolution process represents a novel and effective computational strategy. The use of virtual, rather than real, structured illumination enables the adaptive correction of aberrations, achieving 130 nm lateral resolution at a depth of 180 μm in thick mouse brain tissue. The manuscript is technically sound, with a clearly described algorithmic workflow, comprehensive experimental details, and results that support the applicability of the approach. Overall, this work makes a meaningful contribution to the field and offers a promising computational solution for deep-tissue super-resolution imaging.\nA limitation of the manuscript is the restricted generalizability of the approach. Specifically, the computational adaptive optics framework applies only to virtual structured illumination, which inherently limits imaging speed. Accelerating acquisition would require multifocal scanning strategies, typically involving the capture of hundreds of images, representing a significant weakness compared to hardware AO systems. Additionally, the dual deconvolution algorithm depends on the unique separable structure of the excitation and emission OTFs (Hex and Hem), which is essential for making the problem tractable in this context. Whether this algorithm can be extended to other applications remains unclear and is not discussed in the manuscript.\nAnother limitation lies in the biological demonstrations, which are somewhat limited, featuring only a thin cell and a thick sample with only a few xy planes but lacking any 3D volumetric data. While the lateral resolution enhancement is evident, axial resolution and performance are not discussed. As the deconvolution is currently applied only in 2D, the authors may consider discussing potential extensions of the dual deconvolution algorithm to fully exploit super-resolution capabilities in three dimensions. In addition, Supplementary Figure 3, which presents zebrafish imaging, appears more compelling than Figure 5, as it shows a more pronounced improvement after computational AO correction and effectively demonstrates the method’s ability to address challenging cases with spatially varying aberrations.\nMinor comments:\n1) The data processing speed is not discussed, including the spatially varying iterations for multiple regions.\n2) Line 100, ki is not well defined\n3) Line 104, alpha (numerical aperture?) is not defined.\n4) Fig. 5(g), the dash lines 1-5 seem off from the spine neck.", "reviewer_3": "# Advantages\n1. Interesting idea. The author proposed a dual deconvolution computational aberration correction (AO) algorithm that corrects for aberrations in both excitation and emission in multiphoton deep-tissue structured illumination microscopy (SIM).\n2. No additional hardware required. The proposed method doesn’t require complex wavefront shaping or sensing devices.\n3. Good experimental results. The improvements in aberration correction and resolution are clearly demonstrated.\n4. Code and data are available online, which will benefit the community.\n# Weakness\n1. Computational AO methods generally suffer from low SNR, especially when applied to deep tissue imaging. All experimental results provided by the authors have very high SNR. The paper lacks analysis of how SNR affects the reconstruction quality.\n2. The reconstruction quality in SIM algorithms is highly dependent on the modulation depth or contrast of the illumination pattern, which tends to degrade after propagating through a scattering medium. How do the authors address this issue? And how significantly does this degradation affect the final reconstruction quality?\n3. From a practical standpoint, also as an echo to earlier points, the proposed algorithm is a sophisticated but still blind deconvolution/non-convex optimization method. This kind of approach can be prone to getting trapped in local minima, especially when dealing with complex OTF shapes. The authors highlight that their method requires no prior knowledge of the PSF (as noted in lines 425–426 and 431), and that it uses only simple regularization (line 360, section \"Dual deconvolution algorithm\"). Given this, is there a risk that the algorithm may fail to recover the PSF in certain conditions? I would expect to see examples where the method struggles, such as particular types of excitation/emission PSFs or sample structures that may pose challenges."}, "author_response": "Reviewer’s comment: This article by Choi and Yoon presents a unique approach to achieve super-\nresolution imaging in multiphoton scanning microscopy in aberrating media. The authors use a\nmodified version of a technique known as image scanning microscopy (ISM) that captures images of\nthe emitted point spread function in a laser scanning microscope setup. They demonstrate their novel\napproach of leveraging the diversity of illumination intensity angles and phase shifts present in the\nimaged PSF that naturally emerge during point scanning to perform an alternate version of structured\nillumination. This enables dual-deconvolution computational adaptive optics for the excitation and\nemission light. Their calculated theoretical resolution is 140nm, which is very impressive for distorting\nmedia. In practice, they find closer to 175nm through a scattering phantom and 240nm without\nscattering media. The manuscript is thorough, and presents a compelling case for the novelty and\nimpact of the proposed concept. However, some additional experiments are needed to conclusively\ndemonstrate the technique and to compare with similar current best practices. One notable lack is in\nthe comparison between the described approach and standard ISM. Further discussion of this would\nhelp to strengthen the manuscript. Although the manuscript requires revision, the approach is cutting\nedge and constitutes a significant advance in the use of computational and hardware approaches to\ncreate super-resolution images in living animals. It seems to be of sufficient impact for publication in\nNat Commun and is likely to be of interest to the journal’s readership.\nResponse: We thank the reviewer for the thoughtful and constructive evaluation of our manuscript. We\nare grateful for the recognition of the novelty, rigor, and potential impact of our method, as well as the\nhelpful suggestions for further clarification and validation, particularly the request for a direct\ncomparison with conventional image scanning microscopy (ISM). We performed additional simulations\nand analyses and added a unified theoretical framework for image reconstruction from fluorescence\nresponse function in Supplementary Information (SI). Below, we provide detailed point-by-point\nresponses and indicate how each comment has been addressed. Figures referenced in this response letter\n(Figs. R1-10) are included at the end of this document.\nReviewer’s comment: Given that this paper presents an enhanced version of image scanning\nmicroscopy, a direct comparison in the theoretical and attainable resolution versus typical ISM pixel\nreassignment reconstructions that have been implemented in the past would be relevant to include in\nthe text and figures.\nResponse: We appreciate the reviewer’s insightful suggestion to include a direct theoretical and\nexperimental comparison between the proposed method and standard ISM pixel reassignment. In\nresponse, we have added a comprehensive analysis in SI Section 2, where we present a unified\ntheoretical framework that relates ISM pixel reassignment and SIM spectral synthesis within a common\nrepresentational domain.\nEssentially, the optical transfer functions of ISM and SIM can be expressed as\n𝑘 𝑘\n𝐻 (𝑘) = 𝐻 ( )𝐻 ( ), 𝐻 = 𝐻 ∗𝐻 ,\nISM em ex SIM em ex\n2 2\nwhere * denotes convolution. Therefore, both methods achieve comparable bandwidth extension in\ntheory but differ in their modulation transfer function (MTF) weightings, affecting robustness to noise\nand aberration.\nWe supplemented the theory with simulations and experimental data, as shown in Figures R1–R3:\nFigure R1 illustrates representative image-reconstruction operations of ISM and SIM in both real and\nspatial-frequency spaces, alongside wide-field (WF), laser-scanning (LS), and confocal (CON)\nmodalities for reference.\nFigure R2 provides a direct comparison of five imaging modes (LS, CON, ISM, SIM, AO-SIM) in their\nreconstructed images from the identical fluorescence-response datasets under three conditions: ideal\nsimulation (no aberration), a resolution-target experiment (strong aberration), and a live-cell experiment\n(mild aberration). Under ideal conditions, CON, ISM, and SIM achieve comparable resolution (Fig.\nR2a), whereas their performance diverges in aberrated and noisy environments (Figs. R2b–c). AO-SIM\nimage demonstrates markedly improved image recovery in all cases.\nFigure R3 presents the effective MTFs for ISM and SIM under representative aberrations. Although\nboth cover the same bandwidth, their MTF weightings differ, leading to characteristic attenuation and\ndistortion behaviors. Importantly, the dual-deconvolution algorithm restores the spectrum across the\nentire frequency range, maintaining both resolution and contrast.\nThese results demonstrate that while ISM and SIM offer comparable nominal resolution, our method\nof AO-SIM maintains superior performance in realistic, aberrated environments through explicit\ncorrection of both excitation and emission distortions.\nReviewer’s comment: Lines 122-125— While the details of the dual-deconvolution in the methods\nseems to make sense to this reviewer, the description in the results is harder to follow. It seems\nreasonable that the idea of tuning the phase of each spectrum to enhance the intensity of the synthesized\nspectrum will help with the excitation aberrations. Figure 1 does a decent job of alluding to this\npossibility. However, it is somewhat less clear in the text and figure 1 how the modification of the\nspectral map (as described in lines 137-138) is performed for the emission.\nResponse: We thank the reviewer for pointing this out. We have revised the manuscript to provide the\nintuitive explanation for the emission OTF estimation. The key concept is optical reciprocity. First, we\napply phase-tuning/synthesis to estimate the excitation OTF phase. We then apply the identical\nprocedure to the transposed response matrix, 𝐹T(𝐤 ,𝐤 ) to estimate the emission OTF. Transposing\ni d\nthe response matrix is physically equivalent to interchanging (swapping) illumination and detection\nsides. Therefore, optimizing the synthesized spectrum of 𝐹T allows for the estimation of the emission\nOTF. We have added the following explanation to the main text:\n“An intuitive explanation is to tune the overall phase of each spectrum in Fig. 1d to maximize the total\nintensity of the synthesized spectrum in Fig. 1e, which yields an estimate of 𝐻 . A similar operation\nex\nis conducted on the reciprocal basis (interchanging 𝐤 and 𝐤 ) to estimate 𝐻 . By optical\ni d em\nreciprocity, swapping the illumination and detection ports yields the transpose of the response,\n𝐹T(𝐤 ,𝐤 ), representing the same optics traversed in reverse: the overall phase that previously resided\ni d\non the illumination side now appears on the detection side. Accordingly, we apply the same spectral\nsynthesis and phase tuning to 𝐹T, maximizing the total intensity of the synthesized (transposed)\nspectrum to obtain the estimate of 𝐻 . Iterating these two paths drives the synthesized spectrum\nem\ntoward maximum total intensity (see Methods for detailed mathematical description).”\nReviewer’s comment: Lines 160 and following, 403 and following, and Figure 2— The authors state\nthat “the excitation and emission pupil maps were computationally prepared by the random\nsuperposition of Zernike modes numbering up to 30.” How many random combinations were tested?\nWhat type/depth of tissue aberrations is this hoping to mimic? Were the random aberrations different\nbetween emission and excitation, and how realistic is this combination?\nResponse: We appreciate this important question. A detailed description of the pupil aberration\nmodeling has been added to the SI. We tested the algorithm’s performance over a broad range of random\naberration combinations generated by random superpositions of Zernike modes.\nThe aberrations simulated in Fig. 2 and experimentally demonstrated in Figs. 4f–k are comparable in\nstrength and complexity to those typically encountered in imaging through thinned-skull mouse\npreparations of approximately 30–100 µm thickness, as reported in Nat. Commun. 11, 5721 (2020) and\nNat. Biotechnol. 40, 1663–1671 (2022). For the COS-7 cell and mouse brain tissue experiments shown\nin Fig. 5, the aberrations are relatively mild, consistent with the levels described in PhotoniX 4, 38\n(2023) and Nat. Methods 14, 869 (2017). The zebrafish data (now relocated from the SI to the main\ntext, Fig. 6) exhibit aberration magnitudes in agreement with prior in vivo optical imaging studies (Nat.\nCommun. 10, 3152 (2019)).\nIn general, for one-photon fluorescence imaging, excitation and emission paths experience nearly\nidentical aberrations because of their similar wavelengths. However, in two-photon imaging, these paths\ncan differ substantially. The longer excitation wavelength experiences reduced phase retardation for the\nsame optical-path difference, owing to the inverse wavelength (1/λ) dependence of phase. Therefore,\nexplicit correction of both excitation and emission aberrations is essential to maximize image quality,\nas also emphasized in PhotoniX 4, 38 (2023).\nImportantly, our dual-deconvolution algorithm does not rely on prior knowledge of, or structural\nassumptions about, either PSF. It operates robustly even when excitation and emission aberrations are\ncompletely uncorrelated. To demonstrate this capability, the simulations presented in Fig. 2\nintentionally employed strong, mutually independent aberrations for the excitation and emission pupils.\nThis design choice highlights the algorithm’s ability to recover high-fidelity images under realistic and\nhighly aberrated conditions.\nReviewer’s comment: Are signal loss, noise and tissue scattering (or even aberrations beyond Zernike\nmode 30) likely to play a role in the limitations of the approach? This is mentioned briefly in the\ndiscussion, but additional information could be helpful.\nResponse: Thank you for raising this important point. We have addressed this in SI Sections 4–6, where\nwe systematically evaluate noise sensitivity, aberration strength, and aberration complexity.\n(1) Analysis of SNR and resolution\nWe conducted a SNR analysis while maintaining aberration strength and complexity (Fig. R4).\nIndependent excitation and emission pupil functions were numerically generated using Zernike modes\nup to N=30 (piston and tilt excluded), yielding 505 modes in total. From the corresponding OTFs,\nfluorescence response functions were simulated over a range of noise levels, and AO-2PSIM images\nwere reconstructed via dual deconvolution. As the SNR of the 2PSIM images decreased, the recoverable\nhigh-frequency information in the effective OTFs diminished (as shown in the PTF images in Fig. R4),\nleading to reduced fidelity in both the recovered and corrected PSFs. The algorithm maintains\nreasonable performance down to SNR≈4 (measured from the 2PSIM images), consistent with the\nSNR≈4 estimated for the dataset in Fig. 4h (main text). Below this level, the fidelities of both the\nrecovered target images and the PSFs degrade.\nFourier ring correlation (FRC) curves obtained from 2PSIM images at various SNR conditions show\nresolution degradation with increasing noise before aberration correction (Fig. R5a). In contrast, FRC\ncurves from AO-2PSIM images mostly demonstrate higher cut-off frequencies (Fig. R5b). The direct\nvisualization of these cut-off frequencies (Fig. R5c) shows the bandwidth recovery capability of the\ndual deconvolution method.\nTo further validate this, we performed 10 independent simulation runs at each SNR level and calculated\nSNR-enhancement factors and resolutions. The results show that our method improves SNR by up to\n8-fold (Fig. R5d), and super-resolution (cut-off frequency > 3αkₑₘ, i.e., resolution < λₑₘ/3α) is\nachievable at SNR≈4. These quantitative SNR-resolution correlations clarify the performance of dual\ndeconvolution across diverse noise and aberration conditions.\n(2) Analysis of aberration strength\nTo quantify dual-deconvolution performance as aberration strength increases, we conducted a\nsimulation study under ideal noise-free conditions using Zernike modes up to N=30 (Fig. R6). The\nZernike-coefficient magnitudes were systematically scaled to control aberration strength, and\nreconstruction quality was assessed by the Pearson correlation (PC) between the recovered and ground-\ntruth PSFs. As aberration strength increases, PC values decrease, consistent with reduced high-\nfrequency support in the reconstructed OTFs. At higher aberration levels, the loss of high-frequency\ncontent primarily limits the achievable resolution and renders fine structural details.\n(3) Analysis of aberration complexity\nWe extended the analysis to aberration complexity by increasing the number of Zernike modes from\nN=30 to N=120 (Fig. R7). As complexity increased, the recovered OTFs exhibited increasingly\ncomplex phase structures, and PC between reconstructed and ground-truth images steadily declined,\nindicating a reduction in resolution and accuracy of retrieved OTFs and PSFs.\nReviewer’s comment: Figure 3 and lines 208-219— The bead separation measurement is insightful,\nbut what was the experimentally measured size of the sub-diffraction beads? How close does this come\nto the theoretical predictions?\nIn Fig. 3, we used an excitation wavelength of 𝜆 = 850 nm to image 100-nm-diameter polystyrene\nex\nbeads labeled with Rhodamine B (peak emission 𝜆 = 567 nm). For multiphoton order 𝑛 = 2 and\nem\n1 1\nNA 𝛼 = 1, the cutoff spatial frequency of an ideal OTF is 𝑘 = 2𝜋𝛼( + ), so the ideal\nc\n𝜆em 𝜆ex/𝑛\n𝜋\nbandwidth-limited (Abbe half-pitch) resolution is 𝑟 = = 𝜆 𝜆 /2𝛼(𝑛𝜆 +𝜆 ) ≈ 121 nm.\nideal ex em em ex\n𝑘c\nIn practice, however, the effective passband is set by sampling: the Nyquist cutoff frequency is 𝑘 =\nsc\n𝜋/𝑟 , with camera’s pixel/scan spacing 𝑟 = 131 nm.\nsc sc\nHere we assume that the post-deconvolution OTF is approximately given by a circular flat-top; the\ncorresponding PSF is a jinc-like, PSF(𝑟) ∝ |𝐽 (𝑘𝑟)/(𝑘𝑟)|.\nin the left panel. The orange curve is the ideal PSF profile assuming a perfectly flattened OTF of radius\n𝑘 , 𝐻 (𝑘) = 1 for |𝑘| ≤ 𝜋/𝑟 . The red curve is a model obtained by convolving a 100-nm-\nsc eff sc\ndiameter bead (disk) with a jinc PSF while allowing the band limit to vary as |𝑘| ≤ 𝜋/𝑟 .\nmea\nAlthough\nthe bead size and the post-deconvolution OTF cannot be uniquely separated from a single profile, we\nestimate that the practical band-limited (Abbe) resolution is 𝑟 ≈ 154 nm from the best fit result.\nmea\nReviewer’s comment: It would be interesting to see the calculated excitation and emission Zernike\nmodes for some of the samples displayed. How stable are the optimized Zernike modes for a given set\nof aberrated images?\nResponse: Usually, pupil functions are represented as a combination of Zernike modes. However, in\nour incoherent imaging framework, both excitation and emission OTFs correspond to the\nautocorrelation of the respective pupil functions. Therefore, Zernike coefficients are not directly\naccessible. Instead, the goal of our algorithm is to reconstruct effective intensity PSF and OTF for\naccurate image correction.\nTo assess stability, we examined the reconstructed PSF/OTF across multiple independent aberrated\nimage sets of the same sample. Then, Pearson correlation coefficients were computed between the\nreconstructed and the corresponding ideal, noise-free and aberration-free ground\ntruth images. Our results show that the reconstructions are highly reproducible: both the PSF shape and\nOTF amplitude/phase remain consistent across repeated measurements or simulations with Pearson-\ncorrelation metric.\nAs shown in Fig. R8a, when SNR exceeds ~4, the reconstructed AO-2PSIM images consistently achieve\nover 95% correlation with the ideal images, indicating effective recovery of super-resolution. However,\nbelow this SNR threshold, the correlation drops significantly.\nTo further validate the accuracy of the PSF recovery, we analyzed Pearson correlation between the\nreconstructed and ground-truth excitation and emission PSFs. As shown in Fig. R8b, these correlations\nsteadily improve with increasing SNR and similarly exhibit a clear transition near the super-resolution\nthreshold. This indicates that our reconstruction pipeline can recover both the spatial distribution of the\nexcitation and emission PSFs with high fidelity under low-to-high SNR conditions.\nIn addition, we quantified the enhancement in image quality achieved by our algorithm by computing\nboth the SNR and signal-to-background ratio (SBR) enhancement factors. These were calculated by\ncomparing AO-2PSIM images to standard 2PSIM images. As shown in Fig. R8c, dual deconvolution\nin AO-2PSIM significantly improves image quality, yielding up to 6-fold increases in SNR and up to\n9-fold increases in SBR under high-SNR conditions. These enhancements provide complementary\nevidence of the algorithm's effectiveness in suppressing background noise and enhancing signal contrast.\nAll statistical analyses were based on up to 30 independent simulation runs per condition to ensure\nrobust conclusions.\nIn summary, our extensive simulation results demonstrate that the dual deconvolution AO-2PSIM\nalgorithm reliably recovers super-resolution information under moderate-to-high SNR conditions (SNR\n> 4), with strong correlations to ideal ground-truth images and PSFs. Although we do not explicitly\nextract pupil functions with Zernike modes from the intensity PSFs, the stability of the recovered\nresponse across varying noise and aberration conditions indirectly supports the stability of the\nunderlying phase corrections.\nReviewer’s comment: Figure S3— Since the idea of isoplanatic patch is mentioned, it could support\nthe concept to calculate the isoplanatic patch size for the image shown. This could be accomplished by\napplying the optimized Zernike modes to neighboring locations.\nResponse: To assess spatial PSF variation across a large field of view, we performed patch-wise PSF\nestimation on zebrafish spinal cord imaging data. Figure R9a shows the full field-of-view AO-2PSIM\nimage, identical to Fig. 6 in the revised main text. This image was acquired from the hindbrain spinal\ncord of a 10-dpf zebrafish at a depth of 180 μm. The field of view was divided into a 7×7 grid of\n20 × 20 μm² patches (Fig. R9b), and excitation and emission PSFs were locally estimated from each\npatch using our dual deconvolution algorithm.\nTo quantify PSF similarity across the field of view, we calculated the mean-relative error (MRE)\nbetween each peripheral PSF and the central reference patch (4,4). The MRE metric is calculated as:\n∑ |PSF (𝑖)−PSF (𝑖)|\n𝑖 center peripheral\nMRE(PSF ,PSF ) = .\ncenter peripheral ∑ |PSF (𝑖)|\n𝑖 center\nThe resulting spatial variation map, displayed as 1 – MRE (Fig. R9c), reveals gradual PSF changes\nacross the field of view, with variations increasing toward the edges and reaching approximately 30%\nat the periphery.\nAs an alternative approach to estimate the isoplanatic patch size, we followed the reviewer’s suggestion\n(Fig. R10). From the orange-dotted boxed region in Fig. R9a, we extracted a 10 × 10 μm² central patch\nand designated it as the ground-truth patch (blue dotted box). Neighboring regions within the boxed\narea were then selected at 4 μm intervals. At each location, we independently estimated the excitation\nand emission PSFs. Each of these PSFs was applied to the original ground-truth region to reconstruct\nAO-2PSIM images. As expected, image sharpness degraded progressively as PSFs farther from the\nground-truth location were used.\nTo quantitatively determine the isoplanatic patch size from this analysis, we calculated the MRE curve\nbetween the ground-truth emission PSF and those retrieved from each neighboring location (Fig. R10b).\nThe resulting 1-MRE profile (red curve) was fitted with a smooth exponential function. The 1/e points\nof the fitted curve were used to estimate the isoplanatic patch size, yielding an approximate value of\n14.8 μm (blue dotted line).\nReviewer’s comment: line 312— A little bit more discussion about other approaches to speed imaging\nin ISM and alternate detectors (e.g. SPAD cameras) prior to the concluding paragraph could help to\nstrengthen the case for future impact of the work.\nResponse: We appreciate the reviewer’s suggestion to expand our discussion of approaches for\naccelerating ISM imaging and incorporating alternative detector technologies. While several of these\nstrategies were already mentioned in the concluding paragraph, such as multifocal illumination and\nSPAD array detectors, we have now revised the text to further highlight these points:\n“Future strategies include enhancing image acquisition speed through methods such as parallel imaging\nwith multifocal illumination7,34-36 using a digital micromirror device or a liquid-crystal spatial light\nmodulator19,37-39. Multifocal sparse sampling could further accelerate image acquisition in combination\nwith sparsity-based high-resolution image reconstruction from incomplete measurements40.\nAdditionally, employing SNR-optimized array detectors such as single-photon avalanche diodes\n(SPAD) array detector41 may facilitate capturing raw fluorescence maps at a higher SNR, thereby\nenabling the recovery of the full bandwidth of OTFs.”\nReviewer’s comment: This manuscript presents a solid advancement in multiphoton super-resolution\nimaging through the development of a computational adaptive optics (AO) framework based on virtual\nstructured illumination microscopy (SIM). A central innovation is the dual deconvolution algorithm,\nwhich iteratively corrects excitation and emission aberrations using Wiener-like updates to solve for\nthree unknowns (i.e., object, excitation OTF, and emission OTF), thus decomposing a complex\noptimization problem into well-conditioned, tractable subproblems. While the spectral synthesis\nconcept underlying the two-photon SIM implementation is similar to the recently established k-space\ninterpretation of image scanning microscopy in the single-photon regime, the proposed dual\ndeconvolution process represents a novel and effective computational strategy. The use of virtual,\nrather than real, structured illumination enables the adaptive correction of aberrations, achieving 130\nnm lateral resolution at a depth of 180 μm in thick mouse brain tissue. The manuscript is technically\nsound, with a clearly described algorithmic workflow, comprehensive experimental details, and results\nthat support the applicability of the approach. Overall, this work makes a meaningful contribution to\nthe field and offers a promising computational solution for deep-tissue super-resolution imaging.\nResponse: We greatly appreciate the reviewer’s thoughtful feedback. We are pleased that the reviewer\nrecognizes the novelty and effectiveness of our dual deconvolution algorithm, as well as the significance\nof using virtual structured illumination for adaptive optical aberration correction. We also thank the\nreviewer for acknowledging the clarity of our algorithmic workflow and the technical soundness of our\nexperimental design. These comments reinforce our belief that the proposed method offers a valuable\ncontribution to the field of deep-tissue super-resolution imaging.\nReviewer’s comment: A limitation of the manuscript is the restricted generalizability of the approach.\nSpecifically, the computational adaptive optics framework applies only to virtual structured\nillumination, which inherently limits imaging speed. Accelerating acquisition would require multifocal\nscanning strategies, typically involving the capture of hundreds of images, representing a significant\nweakness compared to hardware AO systems. Additionally, the dual deconvolution algorithm depends\non the unique separable structure of the excitation and emission OTFs (Hex and Hem), which is\nessential for making the problem tractable in this context. Whether this algorithm can be extended to\nother applications remains unclear and is not discussed in the manuscript.\nResponse: We thank the reviewer for this thoughtful and constructive assessment. The comment raises\nthree important aspects: (1) imaging speed, (2) the separable forward model and its generalizability,\nand (3) potential broader applications. We address them in detail below.\n(1) Imaging speed: We agree that imaging speed is a key consideration for the practical use of\ncomputational AO. To provide context, we compared our acquisition rate with that of representative\nhardware AO implementations.\nDirect wavefront-sensing AO typically measures the emission wavefront (e.g., with a Shack–Hartmann\nsensor) using guide stars (well-isolated fluorescent beads or bright point-like emitters), then computes\nDM voltages (often in Zernike modes). In representative biological preparations (e.g., Nat Methods 14,\n869 (2017)), Shack-Hartmann sensor exposures of 0.5–1.5 s per measurement are reported, and the\nmeasure–decompose–apply loop is repeated multiple iterations to reach good correction. The number\nof iterations may change as aberration strength and complexity increases. Typically, end-to-end\ncorrection commonly takes several seconds.\nSensor-less AO, which dispenses with a wavefront sensor, acquires a wide-field image at each iteration\nwhile sweeping many mode coefficients. A commonly cited bound is least 2N+1 wide-field images for\nN Zernike modes, with 7–12 s reported for ~11 modes (Opt. Express 16, 9290 (2008)). In practice, the\nfull optimization can take seconds to minutes depending on the number of modes, signal brightness,\nand metric.\nIn our case, with ~1 ms camera exposure per scan position, a fluorescence response function for\n100×100-pixel FOV (corresponding to ~12 × 12 µm2) is acquired in ~10 s. Moreover, as discussed in\nthe manuscript, multifocal- or instant-SIM schemes (e.g., Nat Methods 9, 749 (2012) and Nat Methods\n10, 1122 (2013)) are directly compatible with our acquisition and can shorten it dramatically; for\nexample, a 10×10 focal-spot array would, in principle, provide a ~100× speed-up at the same exposure\nand SNR.\nConsidering these factors, our computational approach achieves imaging speeds comparable to practical\nhardware AO once measurement overheads are included, while uniquely providing super-resolution\ncapability in an aberrating medium.\nIn fact, the improvement of SNR is a key difference between hardware and computational AO.\nHardware AO physically concentrates light into a tighter PSF and thus improves SNR in the raw image,\nwhile computational AO including our dual-deconvolution method cannot increase the photon budget\nand therefore bounded by the acquisition SNR. We view the two approaches as complementary rather\nthan competing: hardware AO first removes dominant aberrations and raise photon budget, after which\nour spectral synthesis with dual deconvolution recovers residual errors and achieves super-resolution.\nWe detail the limitations of computational AO and emphasize this hybrid perspective in the Discussion:\n“Our proposed computational AO is simpler to implement than hardware AO, but it has drawbacks.\nHardware AO increases raw-image SNR by physically concentrating light into a tighter PSF, whereas\ncomputational AO cannot increase the photon budget; high-spatial-frequency components of the OTF\nare attenuated by aberrations and often buried in noise. Although our dual-deconvolution approach\nrecovers high-frequency content more effectively than conventional deconvolution, its gains are\nultimately bounded by the acquisition SNR. In practice, the two are complementary: hardware AO first\nremoves dominant aberrations and boosts signal, after which spectral synthesis combined with dual\ndeconvolution recovers residual errors and achieves super-resolution.”\n(2) Separable forward model and its generalizability: As the reviewer pointed out, our dual\ndeconvolution leverages a separable forward model:\n𝐹(𝐤 ,𝐤 ) = 𝐻 (𝐤 )𝐻 (𝐤 )𝛤(𝐤),\nd i em d ex i\nwhich ensures the independent estimation of the excitation and emission OTFs, 𝐻 and 𝐻 . This\nex em\nmodel assumes a single object plane at a given depth, with 𝐻 and 𝐻 describing propagation to\nex em\nand from that plane. Multiple scattering and out-of-plane signals are neglected, as is appropriate for\ntwo-photon imaging, which inherently provides optical sectioning.\nThe current formulation considers only main-diagonal elements, corresponding to isoplanatic\nconditions (possessing shift-invariant aberrations). When spatial variability is stronger (anisoplanatic\nfields), 𝐻 and 𝐻 are no longer diagonal matrices. In such cases, a more general three-operator\nem ex\nfactorization is required:\n𝐹(𝐤 ,𝐤 ) = ∑ 𝐻 (𝐤 ,𝐤′)𝛤(𝐤′,𝐤′)𝐻 (𝐤′,𝐤 ).\nd i em d d d i ex i i\n𝐤′,𝐤′\ni d\nTo mitigate this complex factorization problem, we adopted patch-wise processing, with patch sizes\nchosen to match the aberration correlation length. Within each patch, 𝐻 and 𝐻 can be treated as\nex em\nlocally diagonal, allowing stable OTF estimation. This strategy was used in the zebrafish experiment,\nnow presented in Fig. 6 of the revised main text.\n(3) Applicability of the proposed algorithm: We have expanded the Discussion section to clarify\npotential applications beyond the current implementation:\n“Our dual-deconvolution concept is not restricted to SIM. It can support single-molecule localization\nmethods (e.g., STORM/SMLM) and STED microscopy44 by estimating broadened, space-varying PSFs\nfor improved localization. Beyond optics, the dual deconvolution concept extends to an inverse problem\nin systems wherever a linear response function can be measured.”\nReviewer’s comment: Another limitation lies in the biological demonstrations, which are somewhat\nlimited, featuring only a thin cell and a thick sample with only a few xy planes but lacking any 3D\nvolumetric data. While the lateral resolution enhancement is evident, axial resolution and performance\nare not discussed. As the deconvolution is currently applied only in 2D, the authors may consider\ndiscussing potential extensions of the dual deconvolution algorithm to fully exploit super-resolution\ncapabilities in three dimensions. In addition, Supplementary Figure 3, which presents zebrafish\nimaging, appears more compelling than Figure 5, as it shows a more pronounced improvement after\ncomputational AO correction and effectively demonstrates the method’s ability to address challenging\ncases with spatially varying aberrations.\nResponse: We thank the reviewer for raising the perspective of extending the dual deconvolution\nalgorithm to 3D imaging. The simple and straightforward approach would be to perform depth scanning,\nmeasure the fluorescence response function 𝑓(𝐫 ,𝐫,𝑧) at each depth 𝑧, apply dual deconvolution\nd i\nslice by slice, and then construct a 3D stack image. Another possible direction would be to consider a\ntrue 3D dual deconvolution of 𝑓(𝐫 ,𝐫,𝑧) that simultaneously accounts for the 3D OTFs and the\nd i\nvolumetric target distribution. While conceptually simple, this approach would demand substantial\ncomputational and experimental resources, including depth-dependent modeling of both excitation and\nemission PSFs.\nAlternatively, one might attempt to retrieve depth information from a single-depth response function\n𝑓(𝐫 ,𝐫) acquired at the focal plane. This idea is analogous to digital refocusing in coherent imaging\nd i\n(e.g., digital holography). However, in incoherent imaging, the lack of phase information prevents direct\napplication. Moreover, because the 3D intensity PSF is symmetric above and below the focal plane,\nout-of-focus signals cannot easily be distinguished between the +z and –z directions. To break this axial\nsymmetry, engineered PSFs such as those with astigmatic aberration (Proc. Natl. Acad. Sci. U.S.A. 109,\n7175 (2012)) or a double-helix structure (Opt. Lett. 31, 181 (2006)) may be adopted. Nevertheless, the\nimplementation and validation of such 3D extensions are beyond the scope of the present study and\nmerit a dedicated investigation.\nWe have added the following sentence to the Discussion to reflect this future direction:\n“As an intriguing future perspective, the dual deconvolution framework could be extended to three\ndimensions by incorporating depth-dependent PSFs, such as those engineered with astigmatic\naberration42 or a double-helix structure43.”\nIn line with the reviewer’s suggestion, we have also moved the previous Supplementary Fig. S3\n(zebrafish imaging) to the main text as the new Fig. 6.\nReviewer’s comment: (1) The data processing speed is not discussed, including the spatially varying\niterations for multiple regions.\nResponse: In our current implementation, the reconstruction of the large field-of-view zebrafish image\nshown in Fig. 6 (145×145 μm2) took a total of 1,886 seconds in processing 49 patches entirely on the\nCPU. This corresponds to approximately 38 seconds per patch (size of each patch of 20×20 μm2). The\ncomputations were performed on a workstation equipped with an Intel(R) Core(TM) i9-11900K @\n3.50GHz processor and 128 GB RAM. We expect further acceleration through GPU-based\nparallelization and optimized implementation in future work. We added these details to the revised\nmanuscript (in Methods).\nReviewer’s comment: (2) Line 100, ki is not well defined\nResponse: 𝐤 is defined by adding “…with a wavevector 𝐤 ”.\ni i\nReviewer’s comment: (3) Line 104, alpha (numerical aperture?) is not defined.\nResponse: 𝛼 is defined by “…is given by 𝛼𝑘 , where 𝛼 is the numerical aperture and 𝑘 =\nem em\n2𝜋/𝜆 ,…”\nem\nReviewer’s comment: (4) Fig. 5(g), the dash lines 1-5 seem off from the spine neck.\nResponse: The positions of dashed lines are correct. The lower panel shows a saturated image of spine\nneck #1, in which the dimly dendritic spine neck is clearly visualized.\nReviewer’s comment:\n# Advantages\n1. Interesting idea. The author proposed a dual deconvolution computational aberration correction\n(AO) algorithm that corrects for aberrations in both excitation and emission in multiphoton deep-tissue\nstructured illumination microscopy (SIM).\n2. No additional hardware required. The proposed method doesn’t require complex wavefront shaping\nor sensing devices.\n3. Good experimental results. The improvements in aberration correction and resolution are clearly\ndemonstrated.\n4. Code and data are available online, which will benefit the community.\nResponse: We sincerely thank the reviewer for the positive and encouraging feedback. We are pleased\nthat the reviewer recognizes the novelty of our dual-deconvolution computational aberration correction\napproach, as well as its ability to address both excitation and emission aberrations without requiring\nadditional hardware. We also appreciate the acknowledgment of our experimental validation and the\npotential value of openly sharing our code and data. These comments are greatly appreciated and will\nhelp guide future extensions and improvements of our work.\nReviewer’s comment:\n# Weakness\n1. Computational AO methods generally suffer from low SNR, especially when applied to deep tissue\nimaging. All experimental results provided by the authors have very high SNR. The paper lacks\nanalysis of how SNR affects the reconstruction quality.\nResponse: We thank the reviewer for raising this important concern regarding the impact of SNR on\nreconstruction performance, particularly in deep-tissue imaging. To address this point, we performed\nextensive simulations under varying noise conditions and systematically analyzed the robustness of our\nAO-2PSIM algorithm. The results are now included in the Supplementary Information (SI, Sections 4).\nTo investigate dual-deconvolution performance, we conducted simulations across a range of noise\nlevels while maintaining aberration strength and complexity (Fig. R4). Independent excitation and\nemission pupil functions were numerically generated using Zernike modes up to N=30 (piston and tilt\nexcluded), yielding 505 modes in total. From the corresponding OTFs, fluorescence response functions\nwere simulated at various noise levels, and AO-2PSIM images were reconstructed via dual\ndeconvolution. As the SNR of the 2PSIM images decreased, the recoverable high-frequency\ninformation in the effective OTFs diminished (as shown in the PTF images in Fig. R4), leading to\nreduced fidelity in both the recovered and corrected PSFs. The algorithm maintains reasonable\nperformance down to SNR≈4 (measured from the 2PSIM images), consistent with SNR>4 estimated\nfor the dataset in Fig. 4h (main text). Below this level, the fidelities of both the recovered target images\nand the PSFs degrade.\nFourier ring correlation (FRC) curves obtained from 2PSIM images at various SNR conditions show\nresolution degradation with increasing noise before aberration correction (Fig. R5a). In contrast, FRC\ncurves from AO-2PSIM images mostly demonstrate higher cut-off frequencies (Fig. R5b). The direct\nvisualization of these cut-off frequencies (Fig. R5c) shows the bandwidth recovery capability of the\ndual deconvolution method.\nTo further validate this, we performed 10 independent simulation runs at each SNR level and calculated\nSNR-enhancement factors and resolutions. The results show that our method improves SNR by up to\n8-fold (Fig. R5d), and super-resolution (cut-off frequency > 3αkₑₘ, i.e., resolution < λₑₘ/3α) is\nachievable at SNR>4. These quantitative SNR-resolution correlations clarify the performance of dual\ndeconvolution across diverse noise and aberration conditions.\nReviewer’s comment: 2. The reconstruction quality in SIM algorithms is highly dependent on the\nmodulation depth or contrast of the illumination pattern, which tends to degrade after propagating\nthrough a scattering medium. How do the authors address this issue? And how significantly does this\ndegradation affect the final reconstruction quality?\nResponse: Conventional SIM relies on illuminating modulation patterns over a wide field-of-view.\nThese patterns are inherently phase-sensitive, and thus their modulation contrast degrades, and their\nspatial phase is shifted when propagating through deep scattering or aberrating media. Furthermore,\nscattering also increases out-of-focus background signals, which further reduces SNR of raw images.\nHere, the reduction in modulation contrast is described by the modulation transfer function (MTF),\nwhile spatial shifts of the modulation patterns correspond to phase variations in the phase transfer\nfunction (PTF).\nTo mitigate these deteriorating effects, our approach addresses the problem on two complementary\nlevels:\n(1) Physically robust excitation: We employed two-photon excitation with point scanning and\nsynthesized virtual structured illumination, enabling deeper imaging with higher SNR than conventional\none-photon and wide-field illumination SIM. Multiphoton excitation spatially confines emission signal\ninherently to an excitation focal volume such that it effectively suppresses background noise and\nmaintains high-contrast virtual illumination pattern even in scattering samples (Nat. Commun. 16, 5386\n(2025)).\n(2) Computationally robust reconstruction: Even when residual aberrations remain, our dual-\ndeconvolution algorithm explicitly models both excitation and emission aberrations during post-\nprocessing. Conventional SIM deconvolution operates on the synthesized SIM spectrum, 𝐼̃ (𝐤),\nSIM\nwhich suffers from deep local minima in the effective MTF, 𝐻̃ (𝐤), caused by aberrations. It leads\nSIM\nto severe noise amplification during normalization (see Fig. R5). In contrast, our dual-deconvolution\nmethod deconvolves the fluorescence response function, 𝐹(𝐤 ,𝐤 ) before SIM synthesis, using a two-\nd i\ndimensional Wiener-like filter:\n𝐻∗ (𝐤 )𝐻∗ (𝐤 )\nem d ex i\n𝑊 (𝐤 ,𝐤 ) = .\nDCON d i |𝐻 (𝐤 )𝐻 (𝐤 )|2+𝑒\nem d ex i\nThis operation effectively removes phase-induced destructive interference (through conjugate phase\ncorrection) and equalizes amplitude suppression (through MTF compensation). After phase removal,\nthe effective MTF to be corrected reduces to a magnitude-only convolution (see Fig. R5c). Thus, an\neffective MTF for dual deconvolution becomes MTF (𝐤) = (|𝐻 |∗|𝐻 |)(𝐤), thereby avoiding\nDCON em ex\nphase-induced nulls that plague conventional single-OTF-based deconvolution on 𝐼̃ . Consequently,\nISM\nthe dual-deconvolution approach enhances both contrast and stability of the reconstructed SIM images,\neven under significant aberration or scattering.\nReviewer’s comment: 3. From a practical standpoint, also as an echo to earlier points, the proposed\nalgorithm is a sophisticated but still blind deconvolution/non-convex optimization method. This kind of\napproach can be prone to getting trapped in local minima, especially when dealing with complex OTF\nshapes. The authors highlight that their method requires no prior knowledge of the PSF (as noted in\nlines 425–426 and 431), and that it uses only simple regularization (line 360, section \"Dual\ndeconvolution algorithm\"). Given this, is there a risk that the algorithm may fail to recover the PSF in\ncertain conditions? I would expect to see examples where the method struggles, such as particular types\nof excitation/emission PSFs or sample structures that may pose challenges.\nResponse: Conventional blind deconvolution relies on a single measured image whose spectrum can\nbe expressed as the product of the object spectrum and the OTF. The joint estimation of the object and\nOTF spectra is mathematically underdetermined, since their overlapping frequency components result\nin multiple feasible solutions in the absence of prior knowledge or constraints.\nIn contrast, our dual-deconvolution approach utilizes the fluorescence response function, 𝐹(𝐤 ,𝐤 ),\nd i\nwhich captures the object spectra under multiple 𝐤 . The use of multiple measurements removes the\ni\nambiguity between the object and OTF spectra and transforms the problem into a well-determined\nsystem.\nIf 𝐹(𝐤 ,𝐤 ) is represented in in discrete form as an 𝑁×𝑁 fluorescence response matrix, each\nd i\nelement is given by\n𝐹[𝐤 ,𝐤 ] = 𝐻 [𝐤 ]𝛤[𝐤 +𝐤 ]𝐻 [𝐤 ].\nd i em d d i ex i\nHere, the unknowns are the excitation and emission OTFs, 𝐻 [𝐤 ] and 𝐻 [𝐤 ] (each of length 𝑁),\nem d ex i\nand the object spectrum 𝛤 (of length 2𝑁+1), yielding 3𝑁+1 unknowns. Meanwhile, the number\nof equations equals the total number of elements, 𝑁2 , meaning that the system becomes\noverdetermined for sufficiently structured objects.\nHowever, if the object spectrum lacks diversity, the equation becomes linearly dependent, and the\nsolution becomes underdetermined. For example, in the extreme case of a spatially uniform fluorophore\ndistribution whose spectrum contains only DC component, 𝛤 effectively behaves as an identity\nfunction. In this case, the equations become linearly dependent, and the solution becomes\nunderdetermined. In practice, we note that most biological samples possess rich spatial features that\nensure adequate spectral diversity, making the system well-conditioned.\nIn terms of optimization stability, convergence issues arise primarily under low-SNR conditions,\nespecially in high-frequency regions where the OTF magnitude approaches zero. In such cases,\nappropriate regularization (Wiener-like in our implementation) is crucial to stabilize the inversion and\nsuppress noise amplification. With sufficient SNR and object diversity, however, our simulations and\nexperimental results demonstrate stable convergence without the need for any priors or heavy\nregularization.\nIn response to the reviewer’s suggestion, we performed simulations to demonstrate that our dual-\ndeconvolution method can successfully reconstruct both the PSFs, and the object even when the\naberrations are strong and highly complex. The simulations results presented in Figs. R6-7 support the\nconvergence stability of the dual deconvolution algorithm.\n(1) Analysis of aberration strength\nTo quantify dual-deconvolution performance as aberration strength increases, we conducted a\nsimulation study under ideal noise-free conditions using Zernike modes up to N=30 (Fig. R6). The\nZernike-coefficient magnitudes were systematically scaled to control aberration strength, and\nreconstruction quality was assessed by the Pearson correlation (PC) between the recovered and ground-\ntruth PSFs. As aberration strength increases, PC values decrease, consistent with reduced high-\nfrequency support in the reconstructed OTFs. At higher aberration levels, the loss of high-frequency\ncontent primarily limits the achievable resolution and renders fine structural details.\n(2) Analysis of aberration complexity\nWe extended the analysis to aberration complexity by increasing the number of Zernike modes from\nN=30 to N=120 (Fig. R7). As complexity increased, the recovered OTFs exhibited increasingly\ncomplex phase structures, and PC between reconstructed and ground-truth images steadily declined,\nindicating a reduction in resolution and accuracy of retrieved OTFs and PSFs.\nFig. R1 | Comprehensive comparison of image reconstruction methods from the response matrix.\nThis figure illustrates the fundamental differences between various image reconstruction approaches in\nboth spatial (x-domain) and spatial frequency (k-domain) representations. a, Image reconstruction in x-\ndomain. From the measured response function 𝑓(𝑥 ,𝑥), images are reconstructed via simple kernels:\nd i\nwide-field (WF)—sums over all illumination positions, capturing all emitted light regardless of\nexcitation location; laser-scanning (LS)—sums over all detection positions, integrating signal from\npoint illumination; confocal (CON)—extracts main-diagonal slice, combining spatial filtering in both\nexcitation and detection; ISM (pixel reassignment)—sums along anti-diagonals. b, Spectra of the\ncorresponding Fourier domain images in k-space: WF spectrum—slice along detection axis (𝑘 ); LS\nd\nspectrum—slice along illumination axis (𝑘 ); SIM/CON spectrum—anti-diagonal projection; ISM\ni\nspectrum—main-diagonal slice.\nFig. R2 | Comparative image reconstructions from fluorescence response functions. This\nfigure demonstrates practical implementations of five imaging modes on identical datasets a,\nSimulation dataset (aberration-free, noise-free) reconstructed from a response function, showing\ntheoretical limits of each reconstruction method under ideal conditions. b, Experimental dataset\nfrom Fig. 4f-k (main text). c, COS-7 cell specimen from Fig. 5 (main text). Here, five\nreconstruction modes are compared: laser scanning (𝐼 ), confocal (𝐼 ), pixel-reassignment\nLS CON\n(𝐼 ), SIM spectral-synthesis (𝐼 ), and AO-processed AO-SIM (𝐼 ). In b and c, 𝐼\nISM SIM AO−SIM ISM,σ\ndenotes ISM images reconstructed with an optimal Gaussian loose-confocal filter (σ = 0.3 μm).\nAll scale bars indicate 3 μm.\nFig. R3 | Pupil aberrations and effective transfer functions. a, Randomly generated phase\ndistributions for emission and excitation paths (1D slice shown for clarity), representing typical sample-\ninduced aberrations encountered in deep tissue imaging. Here, the numerical aperture 𝛼 = 1 and\n𝑘 = 𝑘 = 𝑘 . b, Single-pass modulation transfer functions (MTFs) correspond to a, |𝐻 (𝑘)|\nmax em ex em\n𝑘 𝑘\nand |𝐻 (𝑘)|. c, Effective MTFs comparison for ISM (𝐻 (𝑘) = 𝐻 ( )𝐻 ( )), SIM (𝐻 =\nex ISM em ex SIM\n2 2\n𝐻 ∗𝐻 ), and dual deconvolution (|𝐻 (𝑘)| = |𝐻 (𝑘)|∗|𝐻 (𝑘)|). The ideal |𝐻 (𝑘)| and\nem ex DCON em ex ISM\nideal |𝐻 (𝑘)| show effective MTFs of ISM and SIM without aberrations, respectively. Aberrations\nSIM\ncreate deeper local suppressions in |𝐻 | than in |𝐻 | , whereas dual deconvolution alleviates\nISM SIM\nphase-induced minimum and yields a broader, better-conditioned passband.\nFig. R4 | Noise analysis demonstrating robustness of AO-2PSIM reconstruction. This\ncomprehensive analysis evaluates algorithm performance across different noise levels while\nmaintaining constant aberration strength and complexity, with rows a-f showing increasing noise\nlevels from SNR 30.4 to 3.9. The first two columns display 2PSIM images before correction and\nAO-2PSIM images after correction with dual deconvolution, respectively, while the third column\nshows estimated excitation and emission PSFs revealing progressive degradation, and the fourth\ncolumn presents recovered phase transfer functions (PTFs) that indicate noise-dependent\nreconstruction fidelity. Scale bars indicate 5 μm in image panels (a, b, d). Scale bars in PTF panels\nindicate emission/excitation cutoff frequencies of αkₑₓ and αkₑₘ, respectively.\nFig. R5 | Quantitative Fourier Ring Correlation (FRC) analysis of resolution and SNR\nenhancement in AO-2PSIM under varying noise conditions. a, FRC curves obtained from\n2PSIM images at various SNR conditions, where the cut-off 1/7 line indicates the cut-off\nfrequencies that define the resolution limit. b, FRC curves obtained from AO-2PSIM demonstrate\nhigher cut-off frequencies compared to those of 2PSIM at various SNR conditions. c, Direct\nvisualization of cut-off frequencies obtained from a and b. d, Combined bar-line plots showing\nSNR enhancement factors and resolution values, with means and standard deviations obtained\nfrom 10 simulation runs at the same SNR but with different aberration and noise fluctuations.\nFig. R6 | Simulation performance of dual-deconvolution algorithm under increasing\naberration strength. Rows a-e show increasing aberration strength. Columns display Zernike\ncoefficients for independently generated excitation and emission pupil functions, the\ncorresponding pupil phase maps and OTF phase maps, AO-2PSIM reconstructions of a Siemens\nstar target, and recovered PSFs alongside ground truth PSFs. RMS WFE: RMS wavefront error\ncalculated from each pupil phase map. P-V: Peak-to-valley calculated from each pupil phase map.\nP-C: Pearson correlation coefficient between recovered and ground-truth PSFs. Scale bars: AO-\n2PSIM, 1.5 μm; PSF, 5 μm. Scale bars in the pupil and OTF maps indicate the excitation/emission\ncutoffs αkₑₓ and αkₑₘ.\nFig. R7 | Simulation performance of dual deconvolution under increasing aberration\ncomplexity. Panels a-e show simulations with different numbers of Zernike modes, using\nindependently generated excitation and emission pupil functions. Across a-e, simulations were\nrepeated with Zernike modes progressively increased from 30 to 120 while the Zernike coefficient\nmagnitudes were kept constant. Zernike modes up to 30 were used in a, and the RMS wavefront\nerror values for the excitation and emission pupil functions are indicated as inset values. The two\ncolumns to the right display the ground-truth excitation phase transfer function (PTF) and\nmodulation transfer functions (MTF) for excitation and emission, respectively. The adjacent\ncolumn shows the reconstructed AO-2PSIM image. Pearson correlation values calculated between\nreconstructed and ground-truth AO-2PSIM images are shown in the figures. Scale bars in the\nexcitation PTF/MTF represent αkₑₓ. Scale bars in the emission PTF/MTF represent αkₑₘ. Scale bar\nin the AO-2PSIM image represents 5 μm.\nFig. R8 | Performance evaluation of the dual deconvolution algorithm under aberration and noise.\na, Pearson correlation between reconstructed and ideal AO-2PSIM images under various SNR\nconditions. Super-resolution reconstruction becomes reliable (≥95% correlation) when SNR exceeds 4.\nb, Correlation of reconstructed excitation/emission PSFs with ground truth, showing consistent\nimprovement beyond the SNR threshold. c, SNR and SBR enhancement factors (AO-2PSIM vs. 2PSIM)\ndemonstrating up to 6× and 9× improvements, respectively, under high-SNR conditions. All metrics\nare averaged over 30 simulation runs.\nFig. R9 | Estimation of spatially varying PSFs over a large field of view from spinal cords of\nzebrafish. a, Identical AO-2PSIM figure from Fig. 6 (main text) is shown. Panel b shows emission\nPSFs reconstructed from panel a. c, Mean-relative error (MRE) values were computed between\nthe PSF obtained from the central patch and peripheral PSFs in panel b, which is represented as a\n1-MRE map. Each individual patch corresponds to the emission PSF shown in b, represented by\ncolumn and row indices.\nFig. R10 | Evaluation of isoplanatic patch size based on MRE metric. a, A ground-truth patch\nwas first defined, and the reconstruction was performed on this region using excitation/emission\nPSFs estimated from neighboring regions. The ground-truth patch (0 μm offset) is represented in\nblue dotted box. PSFs are then estimated from neighboring patches, which are displayed in the\nfirst and second rows, laterally shifted by 4 μm each. All these PSFs were applied to the original\nground-truth region to reconstruct AO-2PSIM images, which are displayed in the third row. Here,\nscalebars represent 2 μm. b, A 1-MRE plot (red fitted curve) shows the evaluation results of spatial\nPSF variance from a. The 1/𝑒 full width was measured to estimate the rough isoplanatic patch\nsize."}, {"round": 2, "reviewer_comments": {"reviewer_1": "No further comments", "reviewer_2": "The authors have adequately addressed my comments. I recommend its publication.", "reviewer_3": "I appreciate the new results provided by the authors. The proposed method is well-demonstrated through comprehensive analysis and I expect it will be applicable to future SIM applications. All my concerns have been addressed and I suggest accept the manuscript."}, "author_response": ""}]

Decision letter