Double-helix optical point spread function enables real-time mesoscopic 3D functional microangiography in the living mouse brain and skull

Domain: neuroscience
Published: 2026-04-13
Editorial outcome: major_revision

Abstract
Quantitative, volumetric imaging of cerebrovascular networks and microcirculation is essential for understanding brain function. However, rapid mesoscopic 3D imaging remains challenging because of fundamental trade-offs between spatiotemporal resolution, field of view, and sensitivity to functional parameters. Here we present a mesoscopic fluorescence imaging platform featuring a double-helix phase mask for real-time, depth-resolved measurements through the intact mouse skull. The compact phase-mask design is compatible with both laser-scanning and widefield microscopy. Using multifocal laser scanning, we demonstrate real-time volumetric in vivo imaging while discriminating calvarial from cerebral vasculature across 6.6×6.6×0.8 mm3 volume. Beyond high-resolution structural imaging, perfusion time-to-peak values are extracted from the laser scanning configuration while accurate flow velocity/direction information is provided via widefield tracking of fluorescently labeled cells. We demonstrate the platform’s capabilities by analyzing brain-layer-specific perfusion dynamics and vascular topology in glioma-bearing mouse brains, offering unprecedented views for probing cerebrovascular alterations in both physiological and pathological contexts.

Introduction
The brain is a highly vascularized organ, with its vasculature forming a complex hierarchical network that spans from large arteries/veins to fine capillaries, delivering essential oxygen and nutrients to support neuronal activity and metabolic demands1,2. This intricate circulatory network plays a vital role in maintaining brain function and homeostasis3. Vascular alterations are hallmarks of many brain disorders, including glioma4,5, stroke6, and Alzheimer’s disease7. Therefore, high-resolution, large-scale volumetric imaging of the cerebral vasculature is critical for advancing our understanding of cerebral physiology8, neurovascular coupling9, and the pathophysiology of multiple neurological disorders10. In addition, growing evidence suggests that calvarial vessels provide critical pathways linking the skull bone marrow and the brain for molecular and cellular trafficking. However, comprehensive visualization of the entire vascular architecture, with sufficient spatial resolution and functional information, remains a persistent technical challenge. Existing imaging modalities are constrained by trade-offs among spatial resolution, field-of-view (FOV), and sensitivity to functional parameters11, limiting their ability to fully capture the structural complexity and flow dynamics of the cerebral vasculature.

Conventional optical imaging techniques based on point-by-point scanning and detection, such as multiphoton microscopy12,13, offer submicron resolution but are restricted to narrow FOVs and lack full-field access to flow velocity/direction measurements. In contrast, localization-based techniques with widefield recordings have demonstrated the capability to visualize capillary-level microcirculation and quantify flow velocity/direction by localizing and tracking intravenously injected microparticles and labeled cells14,15. However, these approaches fail to deliver precise depth information. In a broader view, ultrasound-based16-19 and optoacoustic localization methods20 have enabled 3D visualization of the cerebral vasculature by utilizing matrix, row-column, and spherical transducer arrays. Nevertheless, the cost and complexity of the hardware design, prolonged acquisition times, and reliance on ultrasound coupling hinder broader applicability of these approaches.

Point spread function (PSF) engineering has emerged as a promising strategy to achieve depth-resolved optical imaging in super-resolution microscopy8,21-24, where axial information is encoded into the PSF shape using, e.g., a cylindrical lens or custom phase masks placed in the detection light path. A prominent example is the double-helix PSF (DH-PSF), which encodes depth into the angular rotation of two spatially separated intensity lobes. By precisely measuring the angular orientation of these lobes, the axial position of the emitter can be determined with high precision. This mechanism enables high-precision 3D localization21,25,26, achieving nanometer-scale localization accuracy across a FOV spanning several tens of microns. However, translating this method to mesoscopic imaging is not straightforward due to the low transmission efficiency and prominent field aberration arising from off-axis distortions and imperfections in optical components.

In this work, we built a depth-sensing module based on a customized double-helix (DH) phase mask with an 18 mm aperture size, providing uniform, high-efficiency DH-PSF generation over 800 μm axial range. By directly imposing spatial-domain constraints within an iteration-based phase retrieval framework, the proposed design yields a 3.25-fold improvement in photon efficiency over conventional Laguerre-Gaussian (LG)-based methods27. Additionally, to mitigate field-dependent aberrations under mesoscopic imaging conditions, we further implemented a spatially adaptive calibration strategy, which reduces depth retrieval uncertainty by 8.2-fold. When combined with multifocal laser scanning28,29, our system achieves a lateral resolution of 8.61 μm and an average depth estimation precision of 3.92 μm, supporting volumetric imaging over a 6.6×6.6×0.8 mm3 FOV, with an effective penetration depth of approximately 300 μm through the intact skull. Notably, the gain in depth information is achieved without sacrificing temporal resolution. Leveraging real-time acquisition, we extracted time-to-peak (TTP) values for individual vessels, capturing dynamic perfusion patterns. By utilizing the DH phase mask in a widefield configuration, we reconstructed flow velocity/direction in both transverse and axial planes by tracking fluorescently-labeled red blood cells (RBCs), offering a comprehensive view of structural and functional organization of cortical vasculature.

We extended the application of this technique to glioma-bearing mouse models, enabling volumetric visualization and quantification of tumor-associated vascular abnormalities, including increased vessel tortuosity, disrupted vascular continuity, and altered flow directionality. By enabling depth-resolved imaging over large FOVs with functional readouts, our approach holds promise for advancing fundamental studies on cerebrovascular diseases.

Results
3D fluorescence mesoscopy enabled by an optimized DH-PSF design

The 3D fluorescence mesoscopy platform features a depth-sensing module, which is implemented by inserting a customized DH phase mask between the objective and the tube lens (Fig. 1a). The phase mask (shown in the inset) was positioned at the back focal plane of the objective, encoding the axial position into the angular orientation of the engineered PSF lobes. The compact and modular design enables integration with different types of microscopy systems. Here we considered both laser scanning, as exemplified by multifocal illumination microscopy30,31, and widefield microscopy (Fig. 1a, see Methods for details). In the widefield mode, a collimated laser beam provides uniform epi-illumination, whereas in the laser scanning mode, the beam is diffracted by a beam-splitting grating to generate an array of foci that are rapidly raster scanned across the sample using an acousto-optic deflector (AOD) synchronized with camera acquisition.

To adapt the DH-PSF strategy for mesoscopic imaging, we redesigned DH phase masks following an iterative phase retrieval framework with spatial domain constraints. Specifically, the target amplitude was reshaped into a bimodal Gaussian profile. This optimization resulted in a 3.25-fold increase in photon efficiency (median 0.918 vs. 0.282) compared to the initial design based on LG modes (Fig. 1c). To evaluate the depth sensing performance across the entire FOV, we employed laser scanning at a fixed scan angle and compared representative images with and without the DH phase mask (Fig. 1d). The intensity was quantified at the focal plane by integrating the signal from both lobes of the DH-PSF. Inserting the DH phase mask resulted in a photon efficiency of 54.38%, which remained highly consistent across the FOV (Fig. 1e). The DH-PSF deliberately splits the photon flux into two angularly modulated lobes, reallocating part of the conventional 2D intensity to encode axial information, a trade-off that enables robust depth localization. We next characterized the lateral resolution of the mesoscopic system by imaging DiD-labeled RBCs sparsely distributed on a microscope slide (Fig. 1f), using a 120×120 raster scan with a 3.3 μm step size. The lateral resolution, quantified by the FWHM of Gaussian fits to individual RBC profiles, had a median value of 8.61 μm (Fig. 1h).

To assess depth sensitivity, we imaged a Cy5.5 fluorescent slide while axially scanning across an 800 μm depth range in 25 μm steps using a motorized stage. DH-PSFs from 289 illumination spots were used to fit a linear relationship between the PSF orientation and depth (Fig. 1i). However, due to field-dependent aberrations, noticeable deviations from the fitted linear function were observed especially at the peripheral FOV. These spatial distortions compromise depth retrieval accuracy if uncorrected. To address this, we implemented a field-dependent calibration strategy. The imaging field was partitioned into a 17×17 grid, and for each region, local linear calibration parameters, slope (k) and intercept (b), were independently determined, followed by linear interpolation to form continuous spatial maps (Fig. 1j). Each emitter’s depth was inferred based on its PSF and the local calibration curve corresponding to its lateral position. For validation, we imaged a tilted Cy5.5 fluorescent slide, which presents a constant axial gradient along the vertical axis. The resulting isodepth contours (white lines) in the reconstructed depth map exhibit improved uniformity when using field-dependent calibration curves compared to uncorrected case (Fig. 1k). Quantitatively, the global depth error decreased from 39.17±31.15 μm to 4.79±1.97 μm (Fig. 1l), underscoring the importance of field-dependent correction for accurate depth retrieval in mesoscopic systems. Furthermore, a systematic Monte-Carlo simulation was conducted to assess the robustness and precision of the DH-PSF localization method against noise and photon count (Supplementary Note 1 and Supplementary Fig. 1). In addition, simulations quantifying the impact of phase mask height profile errors confirmed the high tolerance of the DH-PSF phase mask to manufacturing defects (Supplementary Fig. 2).

Furthermore, we benchmarked our optimized DH-PSF design against a previously reported PSF engineering method that employed cylindrical lens to induce astigmatism at the mesoscale29 (Supplementary Fig. 3). Quantitative analysis revealed that the DH-PSF achieved a 14.9% improvement in photon transmission efficiency compared to the cylindrical lens configuration (0.54 ± 0.037 vs 0.47 ± 0.022, mean ± standard deviation, two-sided paired t-test, P = 0.0071, n = 3 pairs). Moreover, the DH-PSF demonstrated enhanced depth estimation precision, as evidenced by a lower standard deviation in the fitted depth values (3.85 μm vs. 5.58 μm; P = 0.026, two-sided paired t-tests). In addition, the off-axis lobe geometry of the DH-PSF inherently reduces central overlap between adjacent emitters, unlike other PSF engineering strategies32,33 with centrally concentrated signals. Leveraging this structural advantage, we introduced an Alternating Direction Method of Multipliers (ADMM)-based reconstruction approach to resolve spatially overlapping PSFs under high emitter densities, as validated by phantom experiments (Supplementary Note 2 and Supplementary Figs. 4–6).

3D transcranial mapping of fluorescent dye perfusion in the murine brain and skull

The skull bone marrow has been increasingly recognized as an active immunological niche, with calvarial vessels functioning as conduits for molecular and cellular exchange between the skull and brain34,35. Discriminating between skull and cortical vasculature is essential for elucidating the spatial organization and connectivity of these vascular networks. Thus, we applied the proposed DH-PSF-based methodology, combined with multifocal illumination, to capture dye-perfusion dynamics in the murine brain transcranially (Fig. 2a, b). Accurate monitoring of perfusion events across depth requires both high axial and temporal resolution. To boost the volumetric imaging speed, we implemented a multi-scan acquisition strategy14, which enables scanning of 2×2 steps within a single camera exposure by leveraging the significantly higher scanning frequency of the AOD relative to the camera frame rate (Fig. 2c). In this configuration, the AOD operated at 4.8 kHz with a base scan grid of 20×20 steps. Under multi-scan mode, this was effectively expanded to a 40×40 sampling grid per volume, enabling an effective volume rate of 3 Hz (Fig. 2d).

A depth-resolved map of cerebral microcirculation in the murine brain acquired transcranially (Fig. 2e) reveals a natural axial gradient across the cortex, caused by the curvature of the skull and brain. Benefiting from the high-speed acquisition, we generated a TTP map, which highlights spatial heterogeneity in vascular perfusion timing (Fig. 2f). By jointly examining the color-coded depth and TTP maps, brain and skull vessels can be readily distinguished based on their relative depths and perfusion delays. Fig. 2g shows a time-lapse image sequence (3.6 s to 7.6 s post-injection) from three representative regions of interest (ROIs) indicated in Fig. 2e, revealing distinct differences in dye arrival times among vascular compartments. Specifically, the fluorescent dye first enters the cerebral microcirculation via intracranial arteries before propagating toward the extracranial calvarial vasculature (Supplementary Movie 1). Within each ROI, a representative vessel pair was selected, one located at a deeper cortical depth (white arrow) and the other at a more superficial position, corresponding to skull vasculature (magenta arrow). Each vessel pair was quantitatively analysed for axial position and perfusion timing (Fig. 2h). Quantitative analysis across n = 18 paired brain-skull vessel segments from 3 mice confirmed consistent differences: brain vessels were 94.40 ± 13.09 μm (mean ± standard error of the mean) deeper than their skull counterparts (two-sided paired t-test, P = 1.3745 × 10-6), and exhibited 4.22 ± 0.53 s smaller TTP values (two-sided paired t-test, P = 3.7764 × 10-7) (Fig. 2i, j). These results demonstrate the system’s capability for unambiguous differentiation of vascular structures located at different depths, facilitating clear separation of the superficial skull vasculature from deeper brain vessels across the entire cortical surface.

3D multiparametric mapping of cerebrovascular dynamics in the murine brain

The DH-PSF sensing module could also be integrated with widefield microscopy to track intravenously injected DiD-labeled RBCs in mice through a cranial window (Fig. 3a). These sparsely labeled RBCs within the vasculature act as point emitters for depth estimation, appearing as characteristic double-lobed DH-PSFs in raw images (Fig. 3b). Each recording lasted 2 minutes. The labeling density was empirically chosen to balance the temporal resolution required for vascular map reconstruction with the need to minimize PSF overlap. Using a lower RBC labeling density would further reduce the probability of PSF overlap but would necessitate proportionally longer recordings to capture sufficient RBC trajectories and flow information for reconstructing a vascular map of comparable quality.

In contrast to the laser scanning mode, where the positions of scanning spots serve as prior information to facilitate the pairing of DH-PSF lobes, the widefield mode necessitated a robust lobe-clustering algorithm. Thus, we first localized the center of lobes in each frame with sub-pixel resolution, followed by the application of a k-means driven spatial clustering approach. In this step, emitters were paired based on spatial proximity, intensity similarity, and geometric constraints (Fig. 3c, see Methods for details). The orientation of DH-PSFs was mapped to depth for each localized RBC to obtain its 3D coordinates, which serves as the input for the 3D tracking algorithm36. These temporally resolved trajectories enabled the reconstruction of volumetric structural map (Fig. 3d), together with their flow directionality (Fig. 3e) and velocity (Fig. 3f, g) by calculating the emitter displacement between the consecutive frames in 3D. Conventional vascular velocity maps, typically derived from planar inter-frame displacements of tracked particles, yield 2D-projected velocity (Vel2D) fields. By incorporating full 3D inter-frame displacements of emitters, we generated both Vel2D (Fig. 3f) and Vel3D (Fig. 3j) maps, considering both transverse and axial velocities. The 3D structural and flow direction maps also facilitate vessel segmentation and vessel type classification. The intensity map of vasculature was first segmented using an Otsu thresholding and morphological closing. Arterial and venous identities were then recognized based on their distinct topological propagation along the vascular graph37,38. This framework yielded branch-averaged maps of flow velocity (Fig. 3h), and vessel type labeling (Fig. 3i) even when pial vessels cross in close proximity. Apparent direct contact between pial arteries and veins appears when one of them dives underneath the other. Comparative analyses across vessel types (n = 10 paired vessels per type) demonstrated significant discrepancies between Vel2D and Vel3D in arteries, arterioles, venules and veins (all P < 0.001), but not in capillaries (P = 0.0527; two-sided paired t-tests; Fig. 3j). The magnitude of the velocity underestimation was pronounced in vessel types such as arterioles and venules, aligning with their greater tendency for steep axial inclinations39,40. The accuracy of axial velocity estimation is tightly linked to axial resolution. At a frame rate of 200 Hz, the axial localization uncertainty of ~4 μm corresponds to an axial velocity uncertainty of ~0.8 mm/s. Reported axial velocities in the murine brain range from 0 to 12 mm/s41,42. For velocities on the higher end of this range, such uncertainty corresponds to ~6.7% relative deviation, and its impact becomes more pronounced for vessels exhibiting lower axial velocities. Consequently, improving axial velocity estimation will fundamentally require enhancing axial localization precision.

Leveraging functional readouts, such as flow velocity and direction, we could further distinguish vertically oriented arterioles and venules spanning across cortical depths (Supplementary Movie 2). This was accomplished by depicting their axial velocity components with a diverging colormap (Fig. 4a), where red and blue denote upward and downward flow, respectively. Two representative penetrating vessels, one arteriole and one venule, marked with black circles in Fig. 4a, were selected for comparative analysis of axial flow characteristics. Time-lapse DH-PSF images of RBCs flowing inside vessels (Fig. 4b, c, Supplementary Fig. 7) exhibited opposite lobe rotations, reflecting their opposing flow directions. These observations were further validated by time-depth plots (Fig. 4d, Supplementary Fig. 7). A magnified ROI is shown in Fig. 4e, emphasizing the anatomical and functional context of the local vasculature. Based on the flow direction map and automated A/V labeling, vascular compartments were labeled and annotated using distinct colors and line styles (red solid line: pial artery; blue solid line: pial vein; red dashed line: penetrating artery; blue dashed line: penetrating vein), aiding in the identification of vessel type and orientation. The corresponding spatiotemporal dynamics of RBC flow within this ROI are visualized in Supplementary Movie 3, which presents the trajectories of individual particles tracked across consecutive frames, highlighting directional flow and vessel-specific transport. Orthogonal y-z slices (Fig. 4f) depict localized axial velocities, demonstrating the system’s ability to resolve both vascular geometry and depth-resolved hemodynamics of cortical microvasculature.

Morphological and functional alterations in tumor-associated vascular remodeling

Access to both morphological and functional readouts make the proposed method a powerful tool for studying vascular alterations during pathological progression, such as during glioblastoma development. Glioblastoma is the most aggressive primary brain tumor43,44, characterized by rapid growth and extensive vascular proliferation. Understanding its vascular remodeling is critical for developing effective therapeutic strategies targeting the tumor microenvironment and mitigating disease progression.

To investigate vascular alterations associated with glioblastoma progression, we used an orthotopic brain tumor model generated via stereotactic injection of U87-MG cells in the right hemisphere of the mouse brain. Prior to functional imaging, high-resolution T1-weighted anatomical reference images were acquired using a 2D fast low angle shot (FLASH) sequence (Supplementary Fig. 8). Cortical blood flow was recorded at 100 Hz under widefield illumination 5 weeks post tumor cell injection, following skull thinning and intravenous injection of fluorescently labeled RBCs. Fig. 5a illustrates the experimental setup and displays a transverse view of the cortex excised after imaging. The ipsilateral side (tumor-bearing hemisphere) exhibited localized hyperfluorescence (red box), arguably indicative of increased vascular permeability and dye extravasation45, whereas the contralateral (non-tumor) hemisphere, lacking such leakage, served as an internal control (blue box). The 3D vascular map of the ipsilateral side revealed a densely interconnected vascular network (Fig. 5b) while the corresponding yz-projection (Fig. 5c) highlights the distribution of penetrating vessels within the tumor area. Structural (Fig. 5d) and velocity (Fig. 5e) maps of the contralateral (healthy) and ipsilateral (tumor) sides demonstrated a dense, disorganized vascular network with irregular branching patterns on the ipsilateral side. To quantitatively assess vascular density, both ipsilateral and contralateral ROIs were subdivided into 10×10 subregions of equal size. The vessel area fraction (VAF) was calculated for each subregion, generating spatial heatmaps (Fig. 5f). To further characterize the distribution of VAF values, probability histograms were plotted for contralateral (Fig. 5g, top) and ipsilateral (Fig. 5g, bottom) ROIs. Statistical comparison of VAF between the ipsilateral and contralateral hemispheres (Fig. 5h) showed a significant increase on the ipsilateral side (Mann-Whitney U test, two-sided; n = 3 mice; N = 300 ROIs per group; U = 36225; P = 3.58×10-5 (approx.)). The median VAF was higher ipsilaterally than contralaterally (0.3347 vs. 0.3173). The observed increase in median VAF is indicative of enhanced angiogenesis and vessel proliferation, likely driven by the tumor's metabolic demands and the resulting hypoxic microenvironment46,47.

One distinct advantage of the tracking-based method over traditional static imaging is its ability to capture RBC trajectories, which is utilized for subsequent analysis of vascular tortuosity and flow directionality. Vascular tortuosity was defined as the ratio of total trajectory length to the straight-line distance between its endpoints. Figs. 6a and 6b show color-coded tortuosity maps in the contralateral and ipsilateral hemispheres, respectively, where ipsilateral side trajectories appeared notably more convoluted48,49. This observation was then confirmed quantitatively by a lower fraction of near-linear paths (tortuosity ≈ 1.0) and a higher incidence in the 1.2-1.6 range (7.7% vs. 1.4%; Figs. 6c, d) in the ipsilateral side. Individual mouse data detailing the shift in the proportion of highly tortuous trajectories is presented in Supplementary Fig. 9. We next assessed flow directionality by computing the orientation of individual trajectories. Direction maps (Figs. 6e, f) revealed highly variable flow directions in tumor vessels, in contrast to the more coherent patterns in the contralateral side. Notably, we observed counter-propagating flow, with instances where RBCs traveled in opposite directions within the same vessel segment. This phenomenon was predominantly observed in narrow, tortuous vessels50,51, and is exemplified in a time-lapse sequence (Fig. 6g) showing two labeled RBCs flowing in opposite directions through a shared vascular path (Supplementary Movie 4). To further investigate flow directionality, we performed angular trajectory analysis on both contralateral and ipsilateral hemispheres. On the contralateral side, three representative vessel segments (Fig. 6h) displayed unidirectional flow, as evidenced by the single-peaked angular distributions in their polar histograms (Fig. 6i). In contrast, in the ipsilateral side, three representative vessel segments (Fig. 6j) exhibited more complex flow patterns. In a middle-sized vessel (segment p4, diameter = 25.2 µm), a single dominant flow direction was observed (Fig. 6k, left), whereas segments p5 and p6 (Fig. 6k, middle and right) showed pronounced bidirectional flow with peaks spaced approximately 180° apart. Such irregular flow behaviors were absent in the contralateral control regions and may reflect impaired perfusion and loss of vascular polarity, known hallmarks of glioma-associated microcirculatory dysfunction51,52. These findings highlight the capability of DH-PSF enabled imaging platform to resolve alterations in vascular phenotypes associated with tumor-driven vascular remodeling.

Methods
DH-PSF based fluorescence mesoscopy

A customized 3D fluorescence mesoscope was developed (Fig. 1a) to support both laser scanning and widefield illumination modes, allowing flexible switching between laser scanning and widefield imaging configurations. Both configurations share the same detection arm that incorporates a DH phase mask positioned at the exit pupil of the objective lens (CLS-SL, EFL = 70 mm, Thorlabs, USA). Fluorescence emission is collected by the objective, passed through a dichroic mirror (F38-663, Semrock, USA) and DH phase mask, and then focused by a tube lens onto a high-speed CMOS camera (1008×1008 pixels, 4.4 kHz at full resolution, pco.dimax S1, PCO AG, Germany). A long-pass emission filter (LP02-671RU-25, Semrock, USA) was placed in front of the camera to block residual excitation light. The magnification ratio of the imaging system was adjusted between 1.50× and 2.86× to accommodate different imaging targets, using two interchangeable tube lenses: TL1 (AF Micro-Nikkor 105 mm, Nikon, Japan) and TL2 (TTL200MP2, Thorlabs, USA).

Widefield illumination was provided by a 640 nm continuous-wave (CW) laser (MDL-HD-640, CNI, China). The laser beam was delivered via a multimode fiber and collimated to ensure uniform illumination on the sample. For brain imaging experiments, the measured illumination power at the sample plane was 139.2 mW. Laser scanning system was implemented using a 660 nm CW laser (gem 660-500 mW, Laser Quantum, USA). The beam was expanded 4× using a relay lens pair (ACN254-050-A and AC254-200-A, Thorlabs, USA), and its polarization angle was adjusted with a half-wave plate (AHWP10M-600, Thorlabs, USA). Beam steering was performed using a two-dimensional acousto-optic deflector (2D AOD, DTSXY-400-532, AA Opto-Electronic, France), which dynamically scanned the beam across the sample plane. The steered beam was directed onto a 17×17 beam-splitting diffractive optical element (DE-R 243, Holoeye Photonics AG, Germany), producing multiple beamlets with an angular separation of 0.27°. This grid density reflects a practical balance between maintaining adequate SNR for each illumination spot and achieving sufficient spatial sampling to model field-dependent aberrations. Synchronization between the AOD and the camera was achieved via external triggering using a data acquisition board (NI PCIe-6536b, National Instruments, USA). These beams were reflected by the dichroic mirror and focused by the objective lens, forming a 2D lattice pattern on the sample. For brain imaging, the total illumination power across the illumination grid was measured at 24.6 mW, corresponding to an average power of 0.085 mW per individual illumination spot.

Design of DH phase masks

To improve energy confinement in the main lobes and enhance the axial fidelity of DH-PSFs, we implemented an iterative phase retrieval framework that integrates a modified amplitude constraint, as illustrated in Supplementary Fig. 10.

We first construct an initial DH-PSF constructed from an equal-weight superposition of LG modes (l,p), with a fixed amplitude (A=1) and phase determined by the LG model. While this initial estimate provides the desired dual-lobe rotation, it exhibits low transmission efficiency and rapid attenuation.

To address these limitations, we defined a double-peaked Gaussian amplitude distribution as the target intensity profile across 39 discrete axial planes. The iterative phase optimization process begins with the initial LG-based DH beam, which is numerically propagated to defocus plane 𝑧𝑘 (for 𝑘 = 1, 2, …, 39) using Fourier transform (𝐹). At each defocus plane 𝑧𝑘, the resulting complex amplitudes 𝑢(𝑧1),…,𝑢(𝑧39) are computed

𝑢(𝑧𝑘) = 𝐹−1[𝐻(𝑧𝑘)⋅𝐹(𝑒𝑖𝜙0)], (1)

where 𝐻 refers to the propagation transfer function, 𝜙0 denotes the initial phase distribution at the plane 𝑧 = 0, and z indicates the propagation distance.

At each plane 𝑧𝑘, a target complex amplitude 𝑢′(𝑧𝑘) is constructed, where the amplitude is defined as the square root of a double-peaked Gaussian function, i.e.,

|𝑢′(𝑧𝑘)| = √𝐺1(𝑥,𝑦;𝜎)+𝐺2(𝑥,𝑦;𝜎). (2)

Each 𝐺𝑖 (for i = 1 and 2) represents a single 2D Gaussian distribution centered at the expected lobe positions as

𝐺𝑖(𝑥,𝑦,𝜎) = 𝑒𝑥𝑝(− (𝑥−𝑥𝑖)2+(𝑦−𝑦𝑖)2 ), (3)

2𝜎2

where(𝑥𝑖,𝑦𝑖) denote the coordinates of the two main lobes in the original complex amplitude, and σ is the standard deviation, optimized to ensure convergence and efficiency.

The amplitudes of 𝑢(𝑧1),…,𝑢(𝑧39) are replaced with 𝑢′(𝑧1),…,𝑢′(𝑧39), while keeping the original phase distribution. An inverse Fourier transform (𝐹−1) is then applied to retrieve the updated phase distribution in the pupil plane as

𝜙𝑛𝑒𝑤 = 𝑎𝑟𝑔(𝐹−1{𝑢′(𝑧)}). (4)

In each iteration, the complex amplitudes obtained across the 39 planes are superposed with equal weights, and the resulting beam is normalized in amplitude before being fed into the next round of optimization. The amplitude correction constraint is enforced by comparing the current amplitude with the predefined double-peak Gaussian function. Through successive iterations, the beam profile progressively converges toward the desired axial and lateral intensity distribution. Upon convergence, the final phase-only mask is extracted and used for fabrication. To accommodate varying imaging depths and optical architectures, we designed three phase masks, each featuring an aperture size of 18 mm. Detailed design specifications are provided in Supplementary Table 2.

Fabrication of DH phase masks

The DH phase masks were fabricated on a fused silica wafer using a combination of multi-level i-line photolithography and reactive-ion etching (RIE). The process utilized a sputtered Cr film (150 nm) as a hard mask. Each lithographic cycle involved hexamethyldisilazane (HMDS) promotion, AZ1505 spin-coating, alignment (1 μm accuracy), and 9 mJ/cm2 UV exposure using an i-line contact aligner. Pattern transfer into the fused silica was achieved through four sequential RIE cycles using CF4 (15 sccm) and O2 (5 sccm) plasma. The cumulative etch depths of 96 nm, 192 nm, 384 nm, and 768 nm yielded the final 16-level diffractive-holographic phase structure, with a maximum height of 1.5374 μm. The step-by-step fabrication procedures are included as Supplementary Table 3.

Calibration of multifocal laser scanning pattern

Prior to image acquisition, the exact positions of the multifocal laser scanning pattern, comprising 289 excitation foci at each scanning angle, were calibrated. This was achieved by temporarily removing the DH phase mask from the detection path and placing a Cy5.5-labeled fluorescent reference slide at the focal plane. For each frame, a 2D Gaussian smoothing filter was applied to the acquired image to reduce noise, and the center of the illumination spots was localized using intensity-weighted centroid localization.

Image reconstruction

An open-source particle tracking software TrackNTrace (TNT, version 1.03)63, originally developed for single-molecule localization microscopy, was used to localize the main lobes of DH-PSF with sub-pixel resolution. Initial emitter detection and position refinement were performed using wavelet-based multi-scale filtering followed by Gaussian fitting. This software outputs a coordinate set of spot candidates at each subframe, which serves as input for subsequent 3D reconstruction.

For datasets acquired under laser scanning mode, the positions of the 289 excitation foci (obtained during the calibration process) were used as prior information. At each illumination spot, a square ROI with a window size of 45 pixels (~330 μm) was defined. Within each ROI, the presence of a dual-lobed DH-PSF was verified based on its expected lobe geometry and signal symmetry. If a valid DH-PSF was found, the central position, intensity, inter-lobe angle, and inferred depth (via angle-to-depth calibration) of emitters were computed. The final volumetric image was rendered by assigning the intensity at the 3D coordinate of localized emitters across all scanning angles.

Unlike the multifocal illumination, widefield illumination provides no prior information about the positions of DH-PSFs. To accurately and efficiently identify dual-lobed DH-PSF emitters, we adopted a clustering-assisted matching strategy based on Lloyd’s k-means algorithm. For each image frame, all refined emitter localizations were partitioned into spatial clusters, with the number of clusters K adaptively set as a function of total emitter count. This pre-clustering step significantly reduced the combinatorial pairing complexity from O(N2) to O(K∙N2), where N denotes the number of emitters in cluster k. By limiting pairing operations to local neighborhoods, this approach greatly improved computational speed while maintaining high matching fidelity. Within each cluster, candidate emitter pairs were screened based on a set of geometric and photometric constraints, including inter-lobe distance, intensity similarity, and angular orientation. Only emitter pairs that satisfied these criteria and were consistent with the expected DH-PSF geometry were retained for depth estimation. This spatially constrained clustering approach preserved reconstruction precision while reducing spurious pairings, and demonstrated improved computational scalability compared to brute-force strategies. The resulting per-frame sparse depth maps were accumulated across time, and a robust per-pixel statistic was applied to suppress outliers and enhance reconstruction fidelity. A pseudocode summary of the main DH-PSF image reconstruction workflow is provided in the Supplementary Information (Supplementary Algorithm 1).

The reconstruction workflow was executed in MATLAB 2023 on a regular computer (Intel i9-13900H CPU, 16.0 GB RAM). Computational times are detailed in Supplementary Table 4.

ADMM-based 3D reconstruction with DH-PSFs

Let 𝐱  ℝ𝑁𝑥×𝑁𝑦×𝑁𝑧 denote the unknown 3D emitter distribution and 𝐲  ℝ𝑁𝑥×𝑁𝑦 the recorded DH-encoded image. The forward model is

𝐲 = 𝐻𝐱+𝒏 = ∑ ℎ𝑧 ∗𝑥𝑧 +𝒏, (5)

where ℎ𝑧 denotes the experimentally calibrated DH-PSF at depth 𝑧, 𝑥𝑧 is the emitter distribution in the corresponding axial plane, * denotes 2D convolution, and n represents measurement noise. To reconstruct 𝐱, we solve the following regularized least-squares problem:

𝑚𝑖𝑛 1 ‖𝐻𝐱−𝐲‖2+𝜆‖𝐱‖1+𝜆𝑇𝑉‖∇𝐱‖1. (6)

𝐱≥0 2

Here, the l1 term promotes sparsity of single-molecule emitters across the 3D volume, while the total variation (TV) term suppresses noise and enforces local spatial consistency. The optimization is carried out with the ADMM64, which alternates between a quadratic update of the volume 𝑥 and soft-thresholding steps for the sparsity and TV terms (Supplementary Note 2). Compared with brute-force or clustering-based approaches, the ADMM formulation offers a key advantage: it can robustly resolve heavily overlapping DH PSFs even under high emitter densities, as validated using phantoms consisting of fluorescent beads at different densities (Supplementary Fig. 5).

Segmentation, vessel graph extraction and A/V labeling

Blood vessels were segmented by first accumulating all tracks persisting for at least three frames, followed by Gaussian blurring and Otsu thresholding of the resulting volume. The segmentation was refined using morphological closing.

A vessel graph was extracted by skeletonizing the segmented volume using the Voreen library’s blood vessel-optimized algorithm, with the bulge-size parameter heuristically optimized to 165. In the resulting graph, nodes represent vessel branching points, and edges correspond to vessel segments between them. Flow direction for each edge was estimated by averaging local flow orientations along the edge’s skeleton voxels, using the local vessel radius as a spatial constraint65.

Automated artery/vein labeling was performed using a recently described global graph optimization method. Briefly, the algorithm leverages the vessel graph topology and known flow directions. Given the high reliability of tracked flow directions, the option to flip flow direction of an edge was disabled. Incorrectly fused A/V crossing points, which can occur when arteries and veins are closely spaced in depth, were resolved by optimizing connectivity to minimize normalized mass flow discrepancies at branch points. Arteries and veins were then classified based on their characteristic diverging (arteries) and converging (veins) flow patterns.

RBC labeling and characterization

RBCs for fluorescence labeling were isolated from donor littermate mice under anesthesia. Fluorescence labeling was performed using a commercial membrane labeling dye called DiD (1,1’-Dioctadecyl-3,3,3’,3’-Tetramethylindodicarbocyanine, D7757, ThermoFisher Scientific, USA), following a previously established protocol66. Briefly, collected blood sample was diluted 10-fold with Dulbecco’s phosphate-buffered saline (DPBS, D8537, Sigma-Aldrich, USA) and centrifuged at 60 ×g for 8 min at room temperature to remove plasma and buffy coat. This washing step was repeated three times to obtain purified RBCs. The resulting RBC pellet was resuspended in 640 µL of DPBS and mixed with 640 µL Diluent C (CGLDIL, Sigma-Aldrich, USA), along with 3.2 µL of DiD dye (50 mg/mL), and incubated for 30 seconds. Unbound dye was removed by four rounds of centrifugation at 60 ×g for 8 min, with DPBS washes between cycles. The DiD-stained RBCs were inspected using a commercial microscope (Cytation C10, BioTek, USA), where the majority of stained RBCs retained their characteristic biconcave disk morphology (Supplementary Fig. 11a). ~2×10⁶ DiD-labeled RBCs were injected into each mouse before imaging. Assuming a total blood volume of 1.5 mL for a 20 g mouse and an average RBC density of ~8.13 × 106 cells/µL67,68, the labeled RBCs account for only ~0.017% of the total RBC population. Additionally, we quantified the number of detectable labeled RBCs over ~25-min observation window (Supplementary Fig. 11b). Linear regression fit revealed a slow decrease trend (clearance rate: -43.41 cells/s), suggesting no acute destabilization or rapid clearance of the labeled RBCs under the experimental conditions.

Animal experiments

All animal experiments were conducted in accordance with the Swiss Federal Act on Animal Protection and approved by the Cantonal Veterinary Office Zurich (ZH182/2023) and by the Animal Experimentation Ethics Committee of Tongji University (TJTJ01925102), Shanghai.

Animal model. C57BL/6J mice (N = 6, 9–14 weeks old, female, Charles River Laboratories, USA) were used to evaluate imaging performance of both systems, with three mice assigned to each configuration. Athymic Foxn1nu mice (N = 3, 13–14 weeks old, female, Charles River Laboratories, USA) were used for tumor studies. Mice were housed in individually ventilated, temperature-controlled cages at 21–23 °C and 40–60% relative humidity under a 12-h light/dark cycle. Pelleted food (3437PXL15, CARGILL) and water were provided ad libitum. Sex was not included as a variable in the experimental design or statistical analysis.

Surgery and anesthesia. Inhalation anesthesia was used with isoflurane (5% for induction, 1.5% for maintenance) with a gas mixture of oxygen (0.1 L/min) and air (0.4 L/min). During surgery and imaging, the mouse was immobilized using a stereotaxic frame (Narishige International Limited, UK), and body temperature was maintained at 37 °C with a feedback-controlled heating system (PhysioSuite, Kent Scientific, USA). Analgesia was provided via subcutaneous injection of buprenorphine (0.1 mg/kg) 30 minutes prior to scalp removal and craniotomy. Laser scanning imaging was performed transcranially after scalp removal, whereas widefield imaging was conducted through a cranial window. This window was created via a 3×3 mm2 craniotomy over the primary somatosensory cortex using a dental drill and sealed with a glass coverslip. Prior to imaging, a catheter filled with phosphate-buffered saline (PBS) was inserted into the tail vein for intravenous injection. Under multifocal illumination configuration, 50 µL of 2 mg/mL fluorescent dye Cy5.5 was administered, whereas for widefield experiments, 100 µL of stained RBC suspension (~2×107 cells/mL) was injected. After the imaging experiment, the mice were euthanized under deep isoflurane anesthesia (5%) for 5 minutes followed by cervical dislocation.

Intracranial tumor implantation. The U87-MG human glioblastoma cell line (Catalog Nr: 300367, CLS Cell Lines Service GmbH, Germany) was cultured in Eagle's Minimum Essential Medium (EMEM, Sigma Aldrich) with 10% FBS (FBS, Sigma Aldrich) in a constant temperature incubator set at 37°C abounding with 5% CO₂. Cells were trypsinized, washed, and resuspended in PBS at 1 × 10⁷ cells/mL for following orthotopic tumor transplantation. Briefly, nude mice were anesthetized with isoflurane (5% for induction and 1.5–2% for maintenance) and received buprenorphine (0.1 mg/kg, s.c.) 30 minutes before surgery. After scalp disinfection, a ~1 cm midline incision was made, and a 0.5-mm burr hole was drilled 1.0 mm posterior to bregma and 1.7–2 mm lateral. A 2 µL suspension (1×10⁶ cells) was injected stereotactically into the brain (0.6–1.2 mm depth) over 15–20 minutes. The hole was sealed with bone wax, and the scalp was sutured. Postoperative care followed institutional guidelines. In this study, orthotopic implantation achieved a 100% survival rate and resulted in tumor growth in all mice. Out of 5 injected mice, we selected 3 mice with tumors growing in the cortex (Supplementary Fig. 8) and excluded the other two because their tumors were situated in the hippocampus. This high level of reproducibility is supported by a standardized cell-culture protocol, stereotactic implantation performed by experienced operators, and precise anatomical targeting guided by a mouse brain atlas. Routine magnetic resonance imaging (MRI) monitoring of tumor development further contributed to the success rate by validating successful engraftment and ensuring that all animals reached an adequate tumor size for imaging analyses. Animals were monitored regularly and euthanized at humane endpoints, including tumor volume ≥1200 mm³ or earlier upon signs of distress (e.g., ulceration, neurological impairment, pain, or >10% body-weight loss). The approved tumor burden limit was not exceeded. Although inter-animal variability in tumor growth rates was observed, longitudinal MRI scans enabled accurate scheduling of downstream experimental procedures at biologically appropriate time points.

Magnetic resonance imaging. Tumor development and anatomical localization were monitored using high-resolution MRI. T1-weighted anatomical reference images were acquired using a two-dimensional fast low angle shot (2D FLASH) sequence. 32 contiguous axial slices were obtained covering the brain from dorsal to ventral. Imaging parameters were as follows: flip angle = 30°, repetition time = 549.7 ms, echo time = 3.5 ms, bandwidth = 44,642.9 Hz, FOV = 15.6×15.218 mm2, matrix size = 312×312, and slice thickness = 0.25 mm, resulting in an effective spatial resolution of 0.050×0.049×0.25 mm3. MRI scans were used to verify tumor engraftment, determine tumor location relative to the imaged FOV in the cortex, and schedule subsequent functional imaging experiments.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents a mesoscopic fluorescence imaging platform based on a double-helix phase mask, enabling real-time depth-resolved measurements of cerebrovascular networks and microcirculation in the living mouse brain through an intact skull. The platform features a compact, generalizable design compatible with both laser scanning and widefield microscopy.\nKey technical advancements include a spatially adaptive calibration strategy to correct field-dependent aberrations, a 3.25-fold improvement in photon efficiency over conventional Laguerre-Gaussian (LG)-based methods via optimized phase mask design, and an ADMM-based reconstruction approach to resolve overlapping point spread functions (PSFs).\nFunctionally, the system achieves a spatial resolution of 8.76 ×8.76 ×3.92 μm3 over a 6.6 ×6.6 ×0.8 mm3 volume (space-bandwidth product of 1.16 ×108) at a 3 Hz volume rate. It extracts hemodynamic parameters: perfusion time-to-peak (TTP) from laser scanning mode and flow velocity/direction from widefield tracking of fluorescently labeled red blood cells (RBCs).\nIn applications, the platform discriminates between calvarial and cerebral vasculature, maps brain-layer-specific perfusion dynamics, and characterizes tumor-associated vascular remodeling (e.g., increased tortuosity, disrupted continuity, bidirectional flow) in glioma-bearing mice, demonstrating utility in both physiological and pathological cerebrovascular studies.\nMajor Revision\nThe manuscript reports a valuable imaging technology with significant advancements in mesoscopic 3D cerebrovascular imaging. Its technical innovations, functional versatility, and in vivo validation in glioma models address critical limitations of existing modalities. However, several substantive issues (detailed below) need clarification and supplement to ensure the rigor, reproducibility, and comprehensiveness of the work before acceptance.\nSpecific Comments and Questions\nTechnical Details and Reproducibility\n1.If this technology is to be widely applied, it is necessary to conduct a thorough assessment of its universality and potential application scope.Can this technology be used to view brain structures other than blood vessels and how to label them?\n2. The phase mask fabrication process mentions using imprint lithography with AZ-1512 photoresist and PDMS replication, but key parameters (e.g., photoresist spin-coating speed, exposure dose, PDMS curing temperature ramp) are missing. Please provide these details to enable reproducibility.\n3. For the spatially adaptive calibration strategy (17×17 grid partitioning), how was the grid size determined? Was a sensitivity analysis conducted to confirm that 17×17 is optimal (e.g., comparing with smaller/larger grid sizes on depth retrieval error)?\n4. The ADMM-based reconstruction approach is introduced to resolve overlapping PSFs, but the manuscript only references Supplementary Note 2 and Fig. 3. Please add a brief description in the main text of the algorithm’s core constraints (e.g., sparsity terms, regularization parameters) and how it outperforms brute-force or other clustering methods.\nFunctional Validation and Quantification\n1. When comparing DH-PSF with cylindrical lens-based PSF engineering (Supplementary Fig. 2), the manuscript reports a 48.3% improvement in photon transmission efficiency and lower depth estimation standard deviation (3.92 μm vs. 5.24 μm). However, it does not specify how many independent experiments were conducted for this comparison. Please provide the sample size (n) and statistical significance tests (e.g., t-test) to support these claims.\n2. In the transcranial dye perfusion experiments (Fig. 2), brain vessels are reported to be 50–150 μm deeper and have 3–6 s earlier TTP than skull vessels. Are these differences consistent across all 6 C57BL/6J mice? Please present pooled data (e.g., box plots) of depth and TTP differences between brain and skull vessels across animals.\n3. For flow velocity measurements in widefield mode (Fig. 3), the manuscript distinguishes between transverse (2D-projected) and total (3D) velocity. What is the average percentage difference between 2D and 3D velocity values across different vessel types (arterioles, capillaries, venules)? This data would highlight the necessity of 3D velocity mapping.\nLimitations and Future Work\n1. The manuscript notes that laser scanning mode underutilizes camera sensor area due to sparse excitation. It proposes wavelength multiplexing or polarization labeling as solutions, but no preliminary data (e.g., feasibility tests, expected improvement in sensor utilization) is provided. Please add a brief discussion of whether these strategies have been tested in pilot experiments.\n2. Deep learning-based reconstruction is suggested as a potential alternative to k-means clustering. Are there any initial results (e.g., comparison of localization accuracy/speed between deep learning and current methods) that can be included to strengthen this proposal?\n3. The system’s effective penetration depth through the intact skull is ~300 μm. How does this compare to other transcranial imaging modalities (e.g., two-photon microscopy with adaptive optics, optoacoustic imaging)? A concise table in the Discussion section comparing penetration depth, FOV, and temporal resolution would contextualize the platform’s performance.\nAnimal Experiments and Ethical Considerations\n1. For the glioma model (U87-MG cell injection), the manuscript states that nude mice were used, but it does not specify the number of mice with successful tumor implantation (e.g., if 3 mice were injected, how many developed tumors suitable for imaging?). Please clarify the tumor implantation success rate and exclude any animals with no/low tumor burden from the analysis.\n2. In the RBC labeling protocol, DiD dye is used at 50 mg/mL, and 3.2 μL is added to the RBC suspension. Is there any evidence of dye toxicity or impact on RBC viability/circulation time (e.g., comparing labeled vs. unlabeled RBC half-life in pilot studies)?\n3. The manuscript mentions that mice were euthanized with ketamine overdose post-imaging. Please confirm that the ketamine dose complies with the Swiss Federal Act on Animal Protection (e.g., dose per gram of body weight) and provide this detail for transparency.\nTypos and Language:\nThe manuscript is well-written, but a few minor corrections are suggested:\nPage 2, Line 74: \"depth-retrival uncertainty\" should be \"depth-retrieval uncertainty\".\nPage 5, Line 257: \"This phenomenon has predominantly manifested\" could be rephrased to \"This phenomenon was predominantly observed\" .", "reviewer_2": "In the manuscript the authors present a novel mesoscopic fluorescence imaging platform which is based on a double-helix phase mask finalized at performing high-speed volumetric imaging of cerebral and calvarial vasculature structures. Double-helix point spread function imaging has been already proposed in microscopy imaging, but here the authors very convincingly extend the technique into the mesoscopic area providing detailed guidance on techniques and strategies for concequent artifacts removal.\nI think that the possibility to perform (as abundantly shown here) in-vivo depth-resolved large-area imaging at the capillary scale through the mouse intact skull is a game-changer of extremely great relevance and offers great potential for probing cerebrovascular alterations as here indicated. It could also find extended use in the cardiovascular field I would imagine.\nThe manuscript is extremely well written. A large amount of carefully given details is present, and it’s very clear in its content, description, and purpose. I found both the exposition and the methodology very clear. Also the language is pretty fluent. The supplementary information is also extremely well done.\nData are very nicely collected and the images appear to be of great quality. I was greatly impressed by the amount of work the authors spent not only in the system characterization but also on how extensively they have tried to demonstrate its use in the biological imaging context.\nAs a minor thing I would explain better the classical double helix point spread function strategy in the introductory section in order to better set a basis for what comes next in the paper. Probably not everybody if familiar with it.\nRegarding the figures they are all excellent in both quality and content description.\nIn Figure 1, in panel “b” it is a little bit difficult to see clearly what is going on.\nPanel “h” indicates a later resolution of 8.76 microns. I would guess this is calculated on one fitted bead. It would be probably better to analyze on more RBCs?", "reviewer_3": "This manuscript presents a new approach to mesoscopic fluorescence imaging by introducing a high-efficiency double-helix point spread functon(DH-PSF) scheme for real-time, depth-resolved microangiography in the intact mouse brain. The work addresses a long-standing challenge in neuro-imaging, namely achieving rapid, volumetric visualization of large-scale cerebrovascular networks without compromising spatial resolution, field-of-view, or functional sensitivity.\nConventional two-photon and localization-based optical methods offer either high spatial resolution or functional flow measurements, but not both simultaneously at mesoscale dimensions. Ultrasound and optoacoustic localization methods provide volumetric access but at the expense of optical resolution, complexity, or transcranial performance. Point spread functon(PSF) engineering has emerged as a promising strategy to achieve depth-resolved optical imaging, where axial information is encoded into the PSF shape, however, particularly DH-PSF has previously not been successfully extended to mesoscopic imaging due to low transmission efficiency and severe field-dependent aberrations.\nThe authors overcome these limitations through a re-designed, 18-mm-aperture, DH phase mask optimised via iterative phase retrieval. Notably, it yields a 3.25-fold gain in photon efficiency. The authors then use a spatially adaptive calibration method to reduce depth-retrieval error by over 8x. When this scheme was combined with multifocal scanning, the system achieved real-time (3 Hz) volumetric imaging over a large 6.6×6.6×0.8 mm³ field at ~9×9×4 µm³ resolution. This enabled reliable discrimination between skull and cortical vessels and quantitative mapping of perfusion dynamics. In widefield mode, the platform further provides 3D flow velocity and direction measurements by tracking labeled RBCs.\nThe ability for both morphological and functional readouts makes the proposed method a powerful tool, which the authors applied to a few compelling use cases. The system was able to visualize depth-resolved perfusion patterns, reconstructed full cortical flow topology, and revealed marked structural and functional abnormalities in glioma vasculature, including increased density, tortuosity, and bidirectional flow.\nOverall, this manuscript introduces a powerful and potentially widely applicable imaging tool which can provide a notable advance for high-throughput mesoscopic neurovascular imaging. As such the manuscript could be a good fit to Nature Communications, however, further clarification and additional validations, as outlined in more detail below, are suggested. Addressing these points would greatly enhance the clarity and impact of the work and bring it closer to the standards expected for a multidisciplinary journal.\nTechnical comments:\n-) I suggest the authors include more systematic quantitative comparisons (resolution, accuracy, imaging depth, speeds, etc) against other state-of-the-art mesoscopic imaging methods, perhaps as a Supplementary Table?\n-) How sensitive is DH-PSF performance to deviations in phase mask height profile or surface roughness? Has this been studied by the authors? Elaborating on this could be useful for further studies and for the general uptake of the technique.\n- ) Was photon efficiency measured across the entire FOV or only on-axis? How does transmission vary off-axis?\n-) How robust is depth estimation to noise and variations in photon counts (signal intensity)?\n-) A quantitative characterization of skull-induced aberrations (e.g., wavefront error maps) would significantly strengthen the claim of robust transcranial performance.\n- ) Regarding the ADMM-based reconstructions and resolution claims: Under high emitter density, how often does the ADMM approach mislocalize overlapped PSFs? How does this perform in more scattering phantoms? Phantom experiment mimicking more realistic conditions of the mouse brain could help to establish realistic error rates.\n-) How sensitive is axial velocity estimation to depth noise, especially given the ~4 μm axial uncertainty?\nBiological demonstrations and findings:\n-) For the skull vs. brain vessel perfusion timing, is the observed separation consistent across all animals? Providing some statistics would strengthen the general findings, also for the other observations presented throughout Figs. 3-6.\n-) Does the overall system performance differ significantly across animals? Has this been tested?\n-) What is the minimal RBC labeling density required for stable flow reconstruction without excessive overlap?\nMinor comments:\n-) Could the authors provide more detail on computational requirements (processing time per volume) so one can assess the scalability for long recordings?"}, "author_response": "Point-by-Point Response to the Reviewers’ Comments\nReviewer #1\nThis manuscript presents a mesoscopic fluorescence imaging platform based on a double-helix\nphase mask, enabling real-time depth-resolved measurements of cerebrovascular networks and\nmicrocirculation in the living mouse brain through an intact skull. The platform features a\ncompact, generalizable design compatible with both laser scanning and widefield\nmicroscopy.\nKey technical advancements include a spatially adaptive calibration strategy to correct field-\ndependent aberrations, a 3.25-fold improvement in photon efficiency over conventional\nLaguerre-Gaussian (LG)-based methods via optimized phase mask design, and an ADMM-\nbased reconstruction approach to resolve overlapping point spread functions (PSFs).\nFunctionally, the system achieves a spatial resolution of 8.76 ×8.76 ×3.92 μm3 over a 6.6 ×6.6\n×0.8 mm3 volume (space-bandwidth product of 1.16 ×108) at a 3 Hz volume rate. It extracts\nhemodynamic parameters: perfusion time-to-peak (TTP) from laser scanning mode and flow\nvelocity/direction from widefield tracking of fluorescently labeled red blood cells (RBCs).\nIn applications, the platform discriminates between calvarial and cerebral vasculature, maps\nbrain-layer-specific perfusion dynamics, and characterizes tumor-associated vascular\nremodeling (e.g., increased tortuosity, disrupted continuity, bidirectional flow) in glioma-\nbearing mice, demonstrating utility in both physiological and pathological cerebrovascular\nstudies.\nMajor Revision\nThe manuscript reports a valuable imaging technology with significant advancements in\nmesoscopic 3D cerebrovascular imaging. Its technical innovations, functional versatility, and\nin vivo validation in glioma models address critical limitations of existing modalities. However,\nseveral substantive issues (detailed below) need clarification and supplement to ensure the\nrigor, reproducibility, and comprehensiveness of the work before acceptance.\nReply: We thank the Reviewer for recognizing the value of our work and insightful suggestions.\nIn the revised manuscript, we have carefully addressed the concerns with additional\nexperiments, simulations, and analyses. Please find below our detailed responses to the specific\ncomments.\nSpecific Comments and Questions\nTechnical Details and Reproducibility\n1. If this technology is to be widely applied, it is necessary to conduct a thorough assessment\nof its universality and potential application scope. Can this technology be used to view\nbrain structures other than blood vessels and how to label them?\nReply: The Reviewer raises an important aspect. While the current study focuses on vascular\nimaging, this choice was mainly motivated by the double-helix (DH) phase mask being\noptimized at a center wavelength of ~700 nm, which facilitates deeper in vivo penetration and\nsignal-to-background ratio (SBR) with far-red fluorescent dyes. Since blood flow naturally\nprovides sparse spatial features and vasculature can be readily labeled with fluorescent dyes,\nwe selected the cerebrovascular network as our initial imaging target. However, the DH-based\nimaging system can potentially be applied to visualize other brain structures, i.e., neurons,\nastrocytes, and glial cells, provided that the phase mask is optimized for a specific excitation\nwavelength and suitable labeling strategies are implemented.\nWith respect to labeling approaches: (1) In widefield imaging scenarios, sparse fluorescent\nlabeling is required to keep individual point emitters distinguishable from neighboring emitters.\nSuch sparsity can be achieved through diluted viral transfections or cell type-specific Cre or\nCreER driver mouse lines [1]. (2) In laser-scanning imaging mode: sparse labeling is not strictly\nrequired, though it might be beneficial for enhancing image contrast. Thus, conventional\ngenetically encoded fluorescent proteins [2] can be used to label cell populations more densely\nthan in widefield mode, which greatly expands the flexibility of the proposed method. However,\nit’s worth mentioning that fluorescent proteins with short emission wavelengths (e.g., ~510 nm\nemission peak for green fluorescent protein) are vulnerable to tissue scattering, whereas longer-\nwavelength variants tend to be dimmer, posing trade-offs in achievable resolution, penetration\ndepth, and signal-to-noise ratio.\nIn the revised manuscript, a discussion on the applicability of the proposed method and the\nlabeling strategies have been added to the Discussion section (page 7, line 333-342)\nReferences:\n[1] Veldman, Matthew B., et al. \"Brainwide genetic sparse cell labeling to illuminate the\nmorphology of neurons and glia with Cre-dependent MORF mice.\" Neuron 108.1 (2020): 111-\n127.\n[2] Chan, Ken Y., et al. \"Engineered AAVs for efficient noninvasive gene delivery to the central\nand peripheral nervous systems.\" Nature neuroscience 20.8 (2017): 1172-1179.\n2. The phase mask fabrication process mentions using imprint lithography with AZ-1512\nphotoresist and PDMS replication, but key parameters (e.g., photoresist spin-coating speed,\nexposure dose, PDMS curing temperature ramp) are missing. Please provide these details to\nenable reproducibility.\nReply: We apologize for incorrect labeling of the fabrication method. The DH phase masks\nwere fabricated using multi-level i-line photolithography and reactive-ion etching (RIE) on a\nfused silica wafer. We have revised the corresponding contents in the Methods section (page 9,\nline 435-442) to explicitly state this process flow and provide step-by-step fabrication\nprocedure and parameters in the new Supplementary Table 4 for reproducibility.\nStep Work Tool Process\nRinse the wafer in Piranha solution (H SO\n2 4\n2 2\ndrying in a wafer drier for 7 min.\nSputtering Deposit 150 nm Cr thin film on the wafer for\n(ESCRD4) 9 min 15 s.\nAdhesion HMDS vapor prime Apply adhesion promotion film for 20 min at\n3\npromotion (YES-310TA) 115 ℃.\nPlace the wafer on the hot plate for 1 min at\n100 ℃.\nAlign the wafer with the mask on the contact\nMask aligner (EVG\n6200)\nwith a dose of 9 mJ/cm2.\nRinse the wafer in Cr Etchant (TechniEtch\nCr01) for 2 min.\nfor 5 min.\nReactive-ion RIE (Plasmalab 100 Etch the wafer with CF (15 sccm) and O (5\n10 4 2\netching - ICP 380) sccm) at 10 ℃.\nRinse the wafer in Cr Etchant (TechniEtch\nCr01) for 5 min.\nSupplementary Table 4: Step-by-step fabrication procedures of DH phase masks.\n3. For the spatially adaptive calibration strategy (17×17 grid partitioning), how was the grid\nsize determined? Was a sensitivity analysis conducted to confirm that 17×17 is optimal (e.g.,\ncomparing with smaller/larger grid sizes on depth retrieval error)?\nReply: Selection of the 17×17 grid is based on (1) availability of beam-splitting gratings on the\nmarket, and (2) a tradeoff between depth fitting precision and SNR. Due to the limited\ncommercially available options, we selected the off-the-shelf 17×17 grating rather than opting\nfor a customized product. To address the second point, we performed an additional simulation,\nas detailed below.\nThe system’s depth retrieval is based on the model depth = (θ-b)/k, where depth is linearly\nlinked to rotation angle of DH-PSF θ. The accuracy of k (slope) and b (intercept) are the critical\ncalibration parameters that directly influence depth retrieval. We simulated the accuracy of k,\nb matrices under four grating configurations with different diffraction orders (N = 5×5, 11×11,\n17×17, 29×29). Assuming a given k, b distribution as the ground truth, each grating corresponds\nto a different sampling density in the original k/b space. The resulting sampling positions on\nthe k and b matrices are shown in the Response Fig. 1a, b. The final k and b maps were then\nobtained by bicubic interpolation based on these sampled points. Increasing N from 5 to 17\nsignificantly reduces the root mean square error (RMSE) for both k and b (Response\nRMSE RMSE\nFig. 1 c, top and middle panels). This confirms that a denser grid is crucial to accurately capture\nthe non-uniform system aberrations and mitigate underfitting. However, increasing N reduces\nthe illumination power deposited into each spot when the total energy is fixed (Response Fig.\n1c bottom), thereby limiting the fluorescence signals in practical measurements. Moreover,\nhigh grid densities may increase the likelihood of PSF overlap. Taken together, the 17×17\nconfiguration reflects a practical balance between maintaining adequate SNR per illumination\nspot and achieving sufficient spatial sampling to model field-dependent aberrations. We have\nadded corresponding descriptions in the Methods section (page 8, line 391-393).\nResponse Figure 1. Simulations on the impact of the number of diffraction order (N) on\ncalibration precision and fluorescence signal. FL: fluorescence.\n4. The ADMM-based reconstruction approach is introduced to resolve overlapping PSFs, but\nthe manuscript only references Supplementary Note 2 and Fig. 3. Please add a brief description\nin the main text of the algorithm’s core constraints (e.g., sparsity terms, regularization\nparameters) and how it outperforms brute-force or other clustering methods.\nReply: As suggested, we have now added a detailed description of the ADMM-based\nreconstruction in the Methods section, detailing the sparsity and regularization terms (page 10,\nline 480-491). Compared with brute-force or clustering-based approaches, the ADMM\nformulation offers a key advantage: it can robustly resolve heavily overlapping DH PSFs even\nunder high emitter densities, as validated in the newly added supplementary Fig. 5 using\nphantoms consisting of fluorescent beads with different densities. By contrast, brute-force or\nclustering methods typically succeed only when the dual lobes of the DH-PSF remain visually\nseparable. The relevant discussion has been added to the Methods section (page 10, line 491-\n494).\nFunctional Validation and Quantification\n1. When comparing DH-PSF with cylindrical lens-based PSF engineering (Supplementary Fig.\n2), the manuscript reports a 48.3% improvement in photon transmission efficiency and lower\ndepth estimation standard deviation (3.92 μm vs. 5.24 μm). However, it does not specify how\nmany independent experiments were conducted for this comparison. Please provide the\nsample size (n) and statistical significance tests (e.g., t-test) to support these claims.\nReply: As suggested, we repeated the comparison between the DH- and cylindrical lens-based\nmethod in N = 3 independent measurements. To accurately capture the characteristic\nmorphology of the PSF under double-helix phase mask and cylindrical lens, we employed a\nglobal analysis based on net integrated intensity, ensuring a robust estimation of photon\nthroughput across all positions. In the revised Supplementary Fig. 3c, g, we have now included\nthe statistical analyses on both metrics (transmission efficiency and standard deviation (SD) of\ndepth estimation) and two-tailed paired t-tests. For photon transmission efficiency, the DH\nmethod yielded an average transmission efficiency of 0.54 compared to 0.47 for the cylindrical\nlens, with a P-value of 0.0071. For the SD of depth estimation, the DH method achieved an\naverage value of 3.85 μm compared to 5.58 μm for the cylindrical lens, with a P-value of 0.0260.\nThe corresponding descriptions have been added to the Results section (page 3, line 140-144)\nand Supplementary Fig. 3.\n2. In the transcranial dye perfusion experiments (Fig. 2), brain vessels are reported to be 50–\n150 μm deeper and have 3–6 s earlier TTP than skull vessels. Are these differences consistent\nacross all 6 C57BL/6J mice? Please present pooled data (e.g., box plots) of depth and TTP\ndifferences between brain and skull vessels across animals.\nReply: The transcranial dye perfusion experiment (Fig. 2) was repeated in N = 3 C57BL/6J\nmice. For statistical analysis, we selected n = 18 brain vessel and n = 18 skull vessels across 3\nmice and presented the pooled depth and time-to-peak (TTP) data as box plots (new Fig. 2i, j).\nStatistical analysis using paired Student’s t-tests confirmed significant disparities in both\nmetrics. Brain vessels were found to be 94.40 ± 13.09 μm (mean ± SEM) deeper than skull\nvessels (P < 0.0001). Furthermore, brain vessels exhibited a significantly earlier TTP by 4.22\n± 0.53 s (P < 0.0001). The corresponding descriptions in the Results section have been updated\n(page 4, line 175-178).\n3. For flow velocity measurements in widefield mode (Fig. 3), the manuscript distinguishes\nbetween transverse (2D-projected) and total (3D) velocity. What is the average percentage\ndifference between 2D and 3D velocity values across different vessel types (arterioles,\ncapillaries, venules)? This data would highlight the necessity of 3D velocity mapping.\nReply: Indeed, this is an important point. We have now quantified the systematic differences\nin absolute values between the transverse velocity (Vel ) and the total velocity (Vel ) across\n2D 3D\narteries, arterioles, capillaries, venules, veins, as shown in the Response Fig. 2a (which has been\nincorporated as the revised Fig. 3j in the manuscript). For analysis, n = 10 vessels per vessel\ntype were selected from 3 mice.\nThe average percentage underestimation of the true velocity (calculated as (Vel - Vel ) /\n3D 2D\nVel × 100%) is presented in Response Fig. 2b. The underestimation is substantial and non-\n3D\ntrivial, ranging from 24.62 ± 11.23 % in capillaries to 53.07 ± 13.91% in venules. The high\nunderestimation in vessels, such as arteries (48.01 ± 10.36%) and venules, is attributed to the\nincreased contribution of the axial velocity component, which is neglected in 2D projected\nmeasurements.\nThe corresponding descriptions have been added to the Results section (page 5, line 208-212)\nand Fig. 3.\nResponse Figure 2. (a) Comparison of absolute Vel (gray) and Vel (blue) velocities, data\n2D 3D\npresented as box plots. (n = 10 per vessel type, ns, not significant, *P < 0.05, **P < 0.01, ***P\n< 0.001). (b) Bar chart of percentage velocity underestimation across different vessel types.\nLimitations and Future Work\n1. The manuscript notes that laser scanning mode underutilizes camera sensor area due to sparse\nexcitation. It proposes wavelength multiplexing or polarization labeling as solutions, but no\npreliminary data (e.g., feasibility tests, expected improvement in sensor utilization) is provided.\nPlease add a brief discussion of whether these strategies have been tested in pilot experiments.\nReply: Wavelength multiplexing represents an important direction for future work. To evaluate\nits feasibility, we have conducted simulations where 488 nm and 660 nm lasers are combined\nfrom the illumination side with lateral displacement induced by diffraction grating. Due to the\ndispersion of beam-splitting grating, these two wavelengths would be projected at different\npositions on the sample plane, causing their corresponding DH-PSFs to occupy different pixel\nregions on the camera sensor (Response Fig. 3). This configuration effectively doubles the\nutilization efficiency of the camera sensor. Additional considerations, such as customized dual-\nband filters and the uniformity of illumination spots produced by the beam-splitting grating,\nwill need to be addressed in practical implementations.\nPolarization-based detection is not directly relevant for increasing the utilization rate of the\ncamera. Although utilizing a polarization camera [1] or pupil-splitting methods [2] have been\nreported to exploit additional camera pixels and extract information such as the orientation of\nfluorescence emission dipole moments in super-resolution microscopy, these approaches are\nless informative in the mesoscopic regime. At this scale, molecular orientations are not fixed\nwithin the relatively large voxel volume, and therefore polarization-based methods may not\nprovide additional information. Thus, we revised the Discussion section to reflect this (page 7,\nline 347-350).\nResponse Figure 3. Simulated spatial intensity distribution under dual-wavelength excitation\n(λ = 488 nm and λ = 660 nm).\n1 2\nReferences:\n[1] Bruggeman, Ezra, et al. \"POLCAM: instant molecular orientation microscopy for the life\nsciences.\" Nature Methods 21.10 (2024): 1873-1883.\n[2] Zhang, Oumeng, et al. \"Six-dimensional single-molecule imaging with isotropic resolution\nusing a multi-view reflector microscope.\" Nature photonics 17.2 (2023): 179-186.\n2. Deep learning-based reconstruction is suggested as a potential alternative to k-means\nclustering. Are there any initial results (e.g., comparison of localization accuracy/speed\nbetween deep learning and current methods) that can be included to strengthen this proposal?\nReply: We appreciate the Reviewer’s suggestion. At the current stage, we have not\nincorporated a deep-learning module into our pipeline nor performed comparisons between\ndeep-learning methods and current k-means-based approach, mainly due to several practical\nchallenges. In particular, we currently lack paired DH raw data and experimentally acquired\nground-truth labels. A potential future direction is to generate synthetic training datasets with\nperfectly paired inputs and ground-truth outputs. However, several challenges must be\naddressed, such as accurately modeling depth- and wavelength-dependent PSF variations,\nasymmetric lobe morphology, and realistic noise levels. Without carefully considering these\nfactors, it is difficult to assemble sufficiently representative training data, which may in turn\ncompromise the generalizability of supervised models. We have added a detailed discussion of\nthese factors in the revised manuscript and outlined our plan to generate realistic DH-PSF\ndatasets to enable a fair evaluation in future work (page 7, line 355-359).\n3. The system’s effective penetration depth through the intact skull is ~300 μm. How does this\ncompare to other transcranial imaging modalities (e.g., two-photon microscopy with adaptive\noptics, optoacoustic imaging)? A concise table in the Discussion section comparing penetration\ndepth, FOV, and temporal resolution would contextualize the platform’s performance.\nReply: We thank the Reviewer for the valuable suggestion. As suggested, we have added\nSupplementary Table 1 directly comparing the key performance parameters of the proposed\nsystem with state-of-the-art transcranial imaging modalities. We have referenced\nSupplementary Table 1 in the Discussion section to emphasize the advantages of our platform\nfor multi-functional neuroimaging (page 7, line 322-324).\nAnimal Experiments and Ethical Considerations\n1. For the glioma model (U87-MG cell injection), the manuscript states that nude mice were\nused, but it does not specify the number of mice with successful tumor implantation (e.g., if 3\nmice were injected, how many developed tumors suitable for imaging?). Please clarify the\ntumor implantation success rate and exclude any animals with no/low tumor burden from the\nanalysis.\nReply: In this study, orthotopic implantation achieved a 100% survival rate and resulted in\ntumor growth in all mice. Animals were selected based on tumor depth (3 of the 5 injected mice\nwere included). We selected the mice with tumors located in the cortex (Supplementary Fig. 8)\nand excluded the other two because their tumors were situated in the hippocampus. This high\nlevel of reproducibility is supported by a standardized cell-culture protocol, stereotactic\nimplantation performed by experienced operators, and precise anatomical targeting guided by\na mouse brain atlas. Routine MRI monitoring of tumor development further contributed to the\nsuccess rate by confirming successful engraftment and ensuring that all animals reached an\nadequate tumor size for imaging analyses. Although inter-animal variability in tumor growth\nrates was observed, longitudinal MRI scans enabled accurate scheduling of downstream\nexperimental procedures at biologically appropriate time points. The corresponding\ndescriptions have been added to the Methods section (page 12, line 564-574).\n2. In the RBC labeling protocol, DiD dye is used at 50 mg/mL, and 3.2 μL is added to the RBC\nsuspension. Is there any evidence of dye toxicity or impact on RBC viability/circulation time\n(e.g., comparing labeled vs. unlabeled RBC half-life in pilot studies)?\nReply: The DiD staining protocol used in our study was adapted and optimized based on\npreviously published work from our group and others [1-2]. To assess the potential influence\nof DiD staining on RBCs, we further examined the stained cells under a commercial microscope.\nThe majority of DiD-stained RBCs retained their characteristic biconcave disk morphology,\nindicating that the staining process does not substantially induce structural damage. A\nrepresentative brightfield image of DiD-stained RBCs has been added as Supplementary Fig.\n11a. It is also important to note that the DiD-stained RBCs constitute only a small fraction\n(~0.017%) of the total circulating RBC population, minimizing any potential systemic toxicity\nor measurable impact on overall blood health or circulation dynamics in the mouse model. The\ncalculation of the ratio of stained RBCs to native RBCs has been added to the Methods section\n(page 11, line 525-528).\nWith respect to RBC viability/circulation time, previous studies report that the half-life of\nnormal RBCs in mice is ~22.6 days [3]. Direct measurement of the in vivo half-life of DiD-\nlabeled RBCs was not feasible in the present study due to the lack of an animal license for\nlongitudinal experiments (required for mice with long-term cranial windows). As an indirect\nassessment, we quantified the number of detectable labeled RBCs over ~25-min observation\nwindow using our imaging system (Supplementary Fig. 11b). A simple linear regression fit\nrevealed a slow decrease trend (clearance rate: -43.41 cells/s), suggesting no acute\ndestabilization or rapid clearance of the labeled RBCs under the experimental conditions.\nMoreover, our functional imaging experiments were completed within a 4.5-min acquisition\nwindow, during which the number of DiD-stained RBCs can be considered stable.\nThe corresponding descriptions have been added to the Methods section (page 11, line 523-\n531).\nSupplementary Figure 11. Assessment of DiD-stained RBC integrity and stability.\nReferences:\n[1] Zhou, Quanyu, et al. \"Cortex-wide transcranial localization microscopy with fluorescently\nlabeled red blood cells.\" Nature Communications 15.1 (2024): 3526.\n[2] Zhang, Zhenkun, et al. \"Imaging volumetric dynamics at high speed in mouse and zebrafish\nbrain with confocal light field microscopy.\" Nature biotechnology 39.1 (2021): 74-83.\n[3] Dholakia, Urshulaa, et al. \"Determination of RBC survival in C57BL/6 and C57BL/6-Tg\n(UBC–GFP) mice.\" Comparative medicine 65.3 (2015): 196-201.\n3. The manuscript mentions that mice were euthanized with ketamine overdose post-imaging.\nPlease confirm that the ketamine dose complies with the Swiss Federal Act on Animal\nProtection (e.g., dose per gram of body weight) and provide this detail for transparency.\nReply: We apologize for the confusion. The mice used in this study were euthanized under\ndeep isoflurane anesthesia (5%) for 5 min followed by cervical dislocation. We confirmed this\nprocedure complies with the Swiss Federal Act on Animal Protection. The corresponding\ndescription has been updated in the Methods section (page 11, line 553-554).\nTypos and Language:\nThe manuscript is well-written, but a few minor corrections are suggested:\nPage 2, Line 74: \"depth-retrival uncertainty\" should be \"depth-retrieval uncertainty\".\nReply: Thank you for pointing out the typo. It has been corrected to \"depth-retrieval\nuncertainty\" in the revised manuscript (page 2, line 76).\nPage 5, Line 257: \"This phenomenon has predominantly manifested\" could be rephrased to\n\"This phenomenon was predominantly observed\".\nReply: The sentence has been revised as suggested (page 6, line 278).\nReviewer #2\nIn the manuscript the authors present a novel mesoscopic fluorescence imaging platform which\nis based on a double-helix phase mask finalized at performing high-speed volumetric imaging\nof cerebral and calvarial vasculature structures. Double-helix point spread function imaging has\nbeen already proposed in microscopy imaging, but here the authors very convincingly extend\nthe technique into the mesoscopic area providing detailed guidance on techniques and strategies\nfor concequent artifacts removal.\nI think that the possibility to perform (as abundantly shown here) in-vivo depth-resolved large-\narea imaging at the capillary scale through the mouse intact skull is a game-changer of\nextremely great relevance and offers great potential for probing cerebrovascular alterations as\nhere indicated. It could also find extended use in the cardiovascular field I would imagine.\nThe manuscript is extremely well written. A large amount of carefully given details is present,\nand it’s very clear in its content, description, and purpose. I found both the exposition and\nthe methodology very clear. Also the language is pretty fluent. The supplementary information is\nalso extremely well done.\nData are very nicely collected and the images appear to be of great quality. I was greatly\nimpressed by the amount of work the authors spent not only in the system characterization but\nalso on how extensively they have tried to demonstrate its use in the biological imaging context.\nReply: We thank the Reviewer for the positive evaluation and constructive suggestions. Please\nfind below our answers to the specific comments.\nAs a minor thing I would explain better the classical double helix point spread function strategy\nin the introductory section in order to better set a basis for what comes next in the paper.\nProbably not everybody if familiar with it.\nReply: We thank the Reviewer for the valuable suggestion. As suggested, we have added a\nbrief description of the original application and basic performance of double-helix point spread\nfunction (DH-PSF) in super-resolution microscopy to the Introduction section (page 2, line 64-\n67).\nRegarding the figures they are all excellent in both quality and content description.\nIn Figure 1, in panel “b” it is a little bit difficult to see clearly what is going on.\nReply: We apologize for the confusion. Fig. 1b compares DH-PSF intensity patterns at different\naxial positions before (top row) and after (bottom row) applying the optimization algorithm. To\nbetter highlight the improvement in the intensity of PSF after optimization, we intentionally\nused different intensity scale ranges for the two rows (0-0.3 before optimization, 0-1.0 after\noptimization). For clarity, we have revised the figure legend to explicitly state the distinct\nintensity scales used for the top and bottom rows of Fig. 1b (page 16, line 754-756).\nPanel “h” indicates a later resolution of 8.76 microns. I would guess this is calculated on one\nfitted bead. It would be probably better to analyze on more RBCs?\nReply: In the original manuscript, the lateral resolution was calculated from a single bead. To\nprovide a more representative assessment, we have randomly selected n = 15 RBCs across the\nfield of view (FOV) for lateral resolution estimation. The analysis yields a lateral resolution of\n8.61 ± 1.41 μm (mean ± SD), as shown in the updated Fig. 1h (page 16, line 761-762).\nReviewer #3\nThis manuscript presents a new approach to mesoscopic fluorescence imaging by introducing\na high-efficiency double-helix point spread function (DH-PSF) scheme for real-time, depth-\nresolved microangiography in the intact mouse brain. The work addresses a long-standing\nchallenge in neuro-imaging, namely achieving rapid, volumetric visualization of large-scale\ncerebrovascular networks without compromising spatial resolution, field-of-view, or functional\nsensitivity.\nConventional two-photon and localization-based optical methods offer either high spatial\nresolution or functional flow measurements, but not both simultaneously at mesoscale\ndimensions. Ultrasound and optoacoustic localization methods provide volumetric access but\nat the expense of optical resolution, complexity, or transcranial performance. Point spread\nfunction (PSF) engineering has emerged as a promising strategy to achieve depth-resolved\noptical imaging, where axial information is encoded into the PSF shape, however, particularly\nDH-PSF has previously not been successfully extended to mesoscopic imaging due to low\ntransmission efficiency and severe field-dependent aberrations.\nThe authors overcome these limitations through a re-designed, 18-mm-aperture, DH phase\nmask optimised via iterative phase retrieval. Notably, it yields a 3.25-fold gain in photon\nefficiency. The authors then use a spatially adaptive calibration method to reduce depth-\nretrieval error by over 8x. When this scheme was combined with multifocal scanning, the\nsystem achieved real-time (3 Hz) volumetric imaging over a large 6.6×6.6×0.8 mm³ field at\n~9×9×4 µm³ resolution. This enabled reliable discrimination between skull and cortical vessels\nand quantitative mapping of perfusion dynamics. In widefield mode, the platform further\nprovides 3D flow velocity and direction measurements by tracking labeled RBCs.\nThe ability for both morphological and functional readouts makes the proposed method a\npowerful tool, which the authors applied to a few compelling use cases. The system was able\nto visualize depth-resolved perfusion patterns, reconstructed full cortical flow topology, and\nrevealed marked structural and functional abnormalities in glioma vasculature, including\nincreased density, tortuosity, and bidirectional flow.\nOverall, this manuscript introduces a powerful and potentially widely applicable imaging tool\nwhich can provide a notable advance for high-throughput mesoscopic neurovascular imaging.\nAs such the manuscript could be a good fit to Nature Communications, however, further\nclarification and additional validations, as outlined in more detail below, are suggested.\nAddressing these points would greatly enhance the clarity and impact of the work and bring it\ncloser to the standards expected for a multidisciplinary journal.\nReply: We sincerely thank the Reviewer for the positive and constructive evaluation of our\nwork. In the revised manuscript, we carefully addressed the concerns from both technical and\nbiological aspects. Specifically, we have provided detailed descriptions of spatial resolution,\nphase mask design, spatially adaptive calibration, and multifocal scanning parameters; a\nsystematic analysis of the impact of photon count and noise levels on depth estimation accuracy;\na comprehensive evaluation of the influence of fluorescent emitter density and scattering\nproperties on the fidelity of ADMM reconstruction, and statistical analyses for biological\nfindings. These clarifications are supported by additional experiments, simulations, and\nanalyses. Please find below our point-by-point responses to the comments.\nTechnical comments:\n-) I suggest the authors include more systematic quantitative comparisons (resolution, accuracy,\nimaging depth, speeds, etc) against other state-of-the-art mesoscopic imaging methods, perhaps\nas a Supplementary Table?\nReply: We fully concur. Accordingly, we have summarized the imaging performance of our\napproach alongside other state-of-the-art mesoscopic imaging modalities (i.e., multi-photon\nmicroscopy, optical coherence tomography, optoacoustic microscopy) in the new\nSupplementary Table 1. We have referenced Supplementary Table 1 in the Discussion section\nto emphasize the unique advantages of our platform for multi-functional neuroimaging (page\n7, line 322-324).\n-) How sensitive is DH-PSF performance to deviations in phase mask height profile or surface\nroughness? Has this been studied by the authors? Elaborating on this could be useful for further\nstudies and for the general uptake of the technique.\nReply: As suggested, we investigated the sensitivity of the double-helix point spread function\n(DH-PSF) to deviations in the phase mask height profile using numerical simulations. Random\nsurface-height perturbations with root-mean-square (RMS) amplitudes of 10 nm, 20 nm, 30 nm,\n40 nm, and 50 nm were introduced to emulate typical fabrication errors. By calculating the\nmaximum intensity of the PSF at various defocus planes and comparing them with the ideal\nerror-free case, we observed negligible performance degradation even at 50 nm height errors\n(RMSE = 3.87%), as shown in the new Supplementary Fig. 2. These results indicate that our\ndesigned DH-PSF exhibits robustness to phase mask manufacturing errors. We have added the\ncorresponding descriptions in the Results section (page 3, line 135-137).\nSupplementary Figure 2. Impact of phase mask manufacturing errors on the maximum intensity\nof DH-PSF across different depths.\n-) Was photon efficiency measured across the entire FOV or only on-axis? How does\ntransmission vary off-axis?\nReply: We thank the Reviewer for raising this point. Photon efficiency was measured across\nthe entire field of view rather than only on-axis by performing a comprehensive, field-wide\nstatistical analysis over the full 17 × 17 spot array (289 spots). To robustly capture the\ncharacteristic morphology of the DH-PSF, we adopted a global metric based on net integrated\nintensity for all spots, which provides a consistent and position-independent estimation of\nphoton throughput.\nAs shown in Response Fig. 4a, the system without phase mask exhibits pronounced spot-to-\nspot intensity fluctuations, with notably reduced intensity in off-axis regions. In contrast, after\nintroducing the DH phase mask, the spot intensity becomes lower but markedly more consistent\nacross the field of view, as reflected by a reduction in the coefficient of variation from 23.05%\nto 18.35%.\nTo further quantify the relative photon transmission across the entire FOV, we computed the\nphoton efficiency on a spot-by-spot basis as the ratio between the integrated intensities with\nand without the DH mask. The resulting spatial efficiency map is shown in Response Fig. 4b.\nWhile slightly higher efficiencies are observed in some peripheral regions, the overall\nefficiency distribution remains highly uniform across the FOV. The global analysis yields a\nmean photon efficiency of 0.544 with a CV of 6.90%, indicating minimal off-axis degradation\nand excellent spatial consistency.\nThe same photon efficiency calculation strategy was consistently applied in Fig. 1e and Fig. S3\n(panel b). The corresponding descriptions have been added to the Results section (page 3, line\n111-112).\nResponse Figure 4. Quantitative evaluation of intensity and photon efficiency. a Spot-by-spot\nintegrated intensity across the full 17 × 17 spot array with (blue) and without (orange) the DH\nphase mask. b Spatial map of photon efficiency, calculated as the ratio of integrated intensity\nwith and without the DH phase mask.\n-) How robust is depth estimation to noise and variations in photon counts (signal intensity)?\nReply: To systematically assess the impact of photon count (N) and noise level (read noise σ\nread\nand background noise B ) on depth estimation precision (σ ), we conducted a three-phase\n0 z\nMonte-Carlo simulation (Supplementary Note 1). We quantified σ across N ranging from 500\nz\nto 10000, σ ranging from 1 to 10 ADU, B ranging from 10 to 500 ADU. The results\nread 0\n(Supplementary Fig. 1) confirm the robustness of DH-PSF depth localization against varying\nnoise characteristics: for photon count (N), precision significantly improved, dropping from an\naverage σ of 9.19 μm (N =500) to 0.93 μm (N =10,000). Regarding camera noise, increasing\nz\nσ from 1 ADU to 10 ADU (at fixed N =2000) resulted in a moderate increase in σ from 2.68\nread z\nμm to 3.88 μm. Increasing B across the full range (10-500 ADU) causes σ to degrade from\n0 z\n2.10 μm to 6.44 μm. Our results confirm that the DH-PSF method maintains high accuracy and\nscales its precision primarily according to theoretical limits, with background noise identified\nas the dominant practical limiting factor. The corresponding contents have been added to\nSupplementary Note 1, Supplementary Fig. 1, and the Results section (page 3, line 133-137).\nSupplementary Figure 1. Systematic analysis of factors influencing DH-PSF depth estimation\nprecision. a-c Dependence of estimation precision (σ) on photon count (N), camera readout\nz\nnoise (σ ), and background noise (B ).\nread 0\n-) A quantitative characterization of skull-induced aberrations (e.g., wavefront error maps)\nwould significantly strengthen the claim of robust transcranial performance.\nReply: As our system does not currently integrate a wavefront-sensing module, we are unable\nto directly measure skull-induced aberrations. Instead, we performed numerical simulations to\nevaluate their influence on DH-PSFs by introducing synthetic aberrations composed of selected\nZernike modes [1-2]. The specific Zernike components and their relative weights (prior to final\nRMS scaling) are summarized in the table below (Response Table 1). Three aberration\nconditions were examined, with total wavefront error RMS values ranging from 2 radians (0.32\nwaves) to 8 radians (1.27 waves).\nCalibration curves relating PSF orientation to depth were simulated using numerical\nconvolution model under these three aberration conditions (Response Fig. 4). The orientation\nangles measured from the aberrated DH- PSF exhibited a minimal deviation from the ideal DH-\nPSF, yielding an RMSE of 0.194°-0.619°. We further performed linear fitting on both the ideal and\nθ\naberrated θ–z curves to estimate the aberration-induced focus shift. The resulting focus shift ranged\nfrom 1.10 to 3.44 μm.\nResponse Table 1: Zernike Mode Components for Aberration Simulation\nNoll Index (j) Zernike Term Relative Amplitude (ΔA, λ)\n10-15 Higher-Order Aberrations 0.01 × randn(6, 1)\nResponse Fig. 4 DHPSF orientation angle error under different aberration severities.\nReferences:\n[1] Yonghyeon Jo et al., Through-skull brain imaging in vivo at visible wavelengths via\ndimensionality reduction adaptive-optical microscopy. Sci. Adv. 8, eabo4366(2022).\n[2] Ji, N., Sato, T. R. & Betzig, E. Characterization and adaptive optical correction of\naberrations during in vivo imaging in the mouse cortex. Proc. Natl. Acad. Sci. U. S. A. 109, 22–\n27 (2012).\n- ) Regarding the ADMM-based reconstructions and resolution claims: Under high emitter\ndensity, how often does the ADMM approach mislocalize overlapped PSFs? How does this\nperform in more scattering phantoms? Phantom experiment mimicking more realistic\nconditions of the mouse brain could help to establish realistic error rates.\nReply: We fully agree with the Reviewer that ADMM-based algorithms may underperform\nunder conditions of high emitter density or strong scattering. To quantitatively assess these\neffects, we performed additional phantom experiments that mimic different emitter densities\nand scattering conditions. Specifically, fluorescent beads with a diameter of ~10 µm were\nimmobilized in 0.65% agar to form a random 3D distribution of emitters, allowing PSF overlap\nto occur more generally, with emitters closely spaced both laterally and axially. The phantom\nwas imaged using a widefield microscope equipped with DH-PSF sensing module. The same\nphantom was also imaged with the same widefield system using a high-NA objective (NA =\n0.25) with z-scanning to render the 3D volume.\nDifferent emitter densities were obtained by adjusting the ratio between fluorescent beads and\nagar, yielding final densities of 7.2×103/ml, 7.2×104/ml, and 3.6×105/ml (Supplementary Fig.\n5). As particle density increased, more overlapping PSFs were observed. ADMM accurately\nrecovered emitter locations at densities of 7.2×103/ml and 7.2×104/ml. At the highest density\n(3.6×105/ml), ADMM still reconstructed all bead positions; however, some incorrect\nlocalizations were observed (5/33 beads), especially in regions of PSF overlap. These errors are\nlikely related to the difficulty in balancing the data-fidelity and TV regularization terms,\ncausing the algorithm to overfit local intensity fluctuations or PSF overlaps.\nTo probe the impact of scattering, we modified the phantom by adding intralipid to the agar at\nfinal concentrations of 0, 0.6%, and 1.2% (Supplementary Fig. 6). In the absence of scattering\n(0%) or under low scattering (0.6%), ADMM reconstruction correctly recovered all bead\nlocalizations. Under increased scattering (1.2%), ADMM still detected all beads, although a\nsmall number of incorrect reconstructions were observed (2/9 beads). The reconstruction errors\nare expected to arise from two aspects: (1) the increasing scattering causes the practical PSF to\ndeviate from the ideal PSF used for reconstruction, and (2) stronger scattering decreases the\nSNR of the acquired images.\nCollectively, these results confirm the method’s ability to resolve spatially and axially\noverlapping emitters over a practical range of emitter densities and scattering properties. The\ncorresponding figures and descriptions have been added in Supplementary Fig. 5,\nSupplementary Fig. 6, and Supplementary Note 3.\n-) How sensitive is axial velocity estimation to depth noise, especially given the ~4 μm axial\nuncertainty?\nReply: Given the 200Hz frame rate in our study, the axial localization uncertainty of ~4 μm\ncorresponds to an axial velocity uncertainty of ~0.8 mm/s. Reported axial velocities in the\nmurine brain range from 0 to 12 mm/s. For velocities on the higher end of this range, such\nuncertainty corresponds to ~6.7% relative deviation, while its impact becomes more\npronounced for vessels exhibiting lower axial velocities. Consequently, improving axial\nvelocity estimation will fundamentally require enhancing axial localization precision. These\ndescriptions have been added to the Results section (page 5, line 212-217).\nReferences:\n[1] Tang, Jianbo, et al. \"Improved color Doppler for cerebral blood flow axial velocity\nimaging.\" IEEE transactions on medical imaging 40.2 (2020): 758-764.\n[2] Miller, D. R., et al. \"Multi-modal microvascular cerebral blood flow velocity mapping\nwith 14T single-vessel MRI and optical microscopy in the mouse brain.\" bioRxiv (2023):\n2023-11.\nBiological demonstrations and findings:\n-) For the skull vs. brain vessel perfusion timing, is the observed separation consistent across\nall animals? Providing some statistics would strengthen the general findings, also for the other\nobservations presented throughout Figs. 3-6.\nReply: We thank the Reviewer for the valuable suggestions regarding the statistical validation\nof our findings. Regarding the skull vs. brain vessel perfusion timing, the observed separation\nwas consistent across all three mice imaged. To further validate this observation, we conducted\na statistical analysis on n = 18 skull vessels and n = 18 brain vessels collected from three mice.\nThe corresponding statistical analysis has been included in the revised Fig. 2j. A paired\nStudent’s t-test confirmed that brain vessels are perfused significantly earlier, exhibiting 4.22\n± 0.53 s lower time-to-peak (TTP) values.\nAs suggested, we have supplemented statistical analyses in Figs. 3-6. Specifically, we have\nadded statistical analyses for both Vel and Vel across different vessel types (revised Fig.\n2D 3D\n3j). We also added a Supplementary Fig. 7 that provides expanded evidence by illustrating the\ndifferentiation of microvessels using axial flow dynamics and depth traces, showcasing\nmultiple representative examples of varying vascular structures. And we have included the\nstatistical results for vessel area fraction (VAF) and trajectory tortuosity (as shown in revised\nFig. 5h and Supplementary Fig. 9).\n-) Does the overall system performance differ significantly across animals? Has this been tested?\nReply: As suggested, we conducted a thorough evaluation of the system’s performance,\nincluding lateral resolution and signal-to-noise ratio (SNR), across different mice and surgical\npreparations, and imaging systems. The lateral resolution was quantified by calculating the full-\nwidth-at-half-maximum (FWHM) of smallest resolvable vessels (n = 6 vessels). The SNR of\nDH-PSFs was calculated by averaging the SNR measured from n = 10 illumination spots (laser\nscanning system) and n = 10 stained RBCs (widefield system). The overall system performance\nwas consistent, as summarized below.\nMouse Surgical Imaging Lateral SNR of DH-\nNumber\nstrain procedure system resolution (μm) PSF\nLaser\nmouse 1 transcranial 10.12 ± 0.35 2.02 ± 1.33\nscanning\nLaser\nmouse 2 transcranial 9.86 ± 0.65 2.38 ± 1.24\nscanning\nLaser\nC57BL/6J mouse 3 transcranial 9.75 ± 0.26 1.81 ± 1.97\nscanning\n(n = 6)\nmouse 4 cranial window Widefield 8.98 ± 0.51 1.37 ± 1.03\nmouse 5 cranial window Widefield 8.86 ± 0.46 1.24 ± 0.76\nmouse 6 cranial window Widefield 9.05 ± 0.24 1.46 ± 0.83\nAthymic mouse 7 skull thinning Widefield 9.12 ± 0.27 1.77 ± 0.54\nFoxn1nu\nmouse 8 skull thinning Widefield 8.92 ± 0.34 1.89 ± 0.79\nmice\n(n = 3) mouse 9 skull thinning Widefield 8.96 ± 0.44 1.82 ± 0.78\n-) What is the minimal RBC labeling density required for stable flow reconstruction without\nexcessive overlap?\nReply: In this study, ~2×10⁶ DiD-labeled RBCs were injected into each mouse before imaging.\nAssuming a total blood volume of 1.5 mL for a 20 g mouse and an average RBC density of\n~8.13 × 106 cells/µL [1-2], the labeled RBCs account for only ~0.017% of the total RBC\npopulation. This labeling density was empirically chosen to balance the temporal resolution\nrequired for vascular map reconstruction with the need to minimize PSF overlap.\nUnder these conditions, each frame contained ~100 detectable DiD-labeled RBCs within a\n672×672 px sensor area, corresponding to an average inter-cell spacing of ~67 px (assuming\napproximately uniform spatial distribution). Given that each DH-PSF spans ~12 px (distance\nbetween the center of dual lobes), the probability of PSF overlap is low. To further mitigate\npotential artifacts from rare overlapping events, only DH-PSFs with high matching scores for\ninterlobe distance, intensity symmetry, or angular-orientation correlation are retained for\ndownstream reconstruction.\nReducing the labeling density would further decrease the probability of PSF overlap but would\nnecessitate proportionally longer recordings to capture sufficient RBC trajectories and flow\ninformation for reconstructing a vascular map of comparable quality.\nThe corresponding descriptions have been updated in the Results section (page 4, line 185-189).\nReferences:\n[1] Leal, Joames KF, Merel JW Adjobo-Hermans, and Giel JCGM Bosman. \"Red blood cell\nhomeostasis: mechanisms and effects of microvesicle generation in health and disease.\"\nFrontiers in physiology 9 (2018): 703.\n[2] Huisjes, Rick, et al. \"Squeezing for life-properties of red blood cell deformability.\"\nFrontiers in physiology 9 (2018): 656.\nMinor comments:\n-) Could the authors provide more detail on computational requirements (processing time per\nvolume) so one can assess the scalability for long recordings?\nReply: As suggested, we have added Supplementary Table 5, summarizing the processing time\nfor each step of image reconstruction workflow in both the laser scanning and widefield\nmicroscope systems. We have added the corresponding descriptions in the Methods section\n(page 10, line 478-479)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed most of my previous concerns through detailed explanations and substantial revisions to the manuscript. The additional data provided in the revised paper significantly strengthens the paper’s conclusions.\nHowever, a few minor points remain that require further clarification or technical refinement before the manuscript can be recommended for publication.\nData and Code Availability (Crucial for Reproducibility):\nThe manuscript relies heavily on custom algorithms (ADMM, K-means clustering) and a specialized phase plate design. To ensure the reproducibility of the findings and the utility of this platform for the broader imaging community, the authors should provide more than a \"upon request\" statement.\nI suggest the authors deposit the core reconstruction scripts and a representative raw data sample (e.g., a short image sequence of the double-helix lobes) to a public repository like GitHub or Zenodo. At minimum, a detailed pseudo-code or a flowchart of the ADMM optimization process should be included in the Supplementary Information.", "reviewer_2": "The authors replie to all the comments I had. I also found the replies to the other reviewers' comments very insightful.", "reviewer_3": "I have read the revised manuscript by the authors and found it much improved. My comments and questions were thoroughly answered and addressed. In particular, the more quantitative comparison as well as characterisation of the system performance in the presence of noise add valuable details. The additional experiments statistics for the animal experiments also strengthen the findings of the manuscript. I therefore recommend publication."}, "author_response": "Point-by-Point Response to the Reviewers’ Comments\nReviewer #1\nThe authors have addressed most of my previous concerns through detailed explanations and\nsubstantial revisions to the manuscript. The additional data provided in the revised paper\nsignificantly strengthens the paper’s conclusions.\nHowever, a few minor points remain that require further clarification or technical refinement\nbefore the manuscript can be recommended for publication.\nData and Code Availability (Crucial for Reproducibility):\nThe manuscript relies heavily on custom algorithms (ADMM, K-means clustering) and a\nspecialized phase plate design. To ensure the reproducibility of the findings and the utility of\nthis platform for the broader imaging community, the authors should provide more than a \"upon\nrequest\" statement.\nI suggest the authors deposit the core reconstruction scripts and a representative raw data\nsample (e.g., a short image sequence of the double-helix lobes) to a public repository like\nGitHub or Zenodo. At minimum, a detailed pseudo-code or a flowchart of the ADMM\noptimization process should be included in the Supplementary Information.\nReply: We thank the reviewer for this valuable suggestion regarding reproducibility and\naccessibility of the reconstruction workflow. In response, we have added a detailed pseudo-\ncode for the double-helix image reconstruction pipeline in the Supplementary Information\n(Supplementary Algorithm 1). The ADMM implementation used in this study builds on the\nopen-source DiffuserCam framework; we now cite the corresponding reference in the Code\navailability section (line 584, page 12) and include the GitHub repository link in the Reporting\nSummary. In addition, we have deposited representative raw data collected under both laser-\nscanning and wide-field illumination modes in Zenodo at:\nhttps://doi.org/10.5281/zenodo.18876905. The Data Availability section in the manuscript and\nReporting Summary has been updated accordingly (line 579–581, page 12).\nReviewer #2\nThe authors replie to all the comments I had. I also found the replies to the other reviewers'\ncomments very insightful.\nReply: We thank the Reviewer for the positive feedback on our revision and are sincerely\ngrateful for the constructive suggestions, which have been very helpful in improving the quality\nof this work.\nReviewer #3\nI have read the revised manuscript by the authors and found it much improved. My comments\nand questions were thoroughly answered and addressed. In particular, the more quantitative\ncomparison as well as characterisation of the system performance in the presence of noise add\nvaluable details. The additional experiments statistics for the animal experiments also\nstrengthen the findings of the manuscript. I therefore recommend publication.\nReply: We thank the Reviewer for the positive feedback on our revision. We are grateful for\nthe valuable suggestions in the last round, which have been very helpful in improving the\nquality of this work."}]

Decision letter