Fluorescence microscopy (FM) can provide both in vitro and in vivo imaging of biological tissues as well as their functional dynamics, which is observable with high spatial and temporal resolution [1–3]. Moreover, benefiting from the various choices of fluorescent indicators for neuronal activities, FM has been widely used for in vivo brain imaging [4–9]. However, because a large number of neurons are spread in the brain region [10], brain imaging must record neuronal activities with a large field of view (FOV) simultaneously [11–13]. For example, three-dimensional imaging of the zebrafish brain should cover a volume of 800μm×600μm×300μmwith subneuron resolution.

Different from confocal [14], spinning disk confocal [15], multiphoton [16], and light-field [17] FMs, selective plane illumination microscopy [18–20] (also called light-sheet microscopy) has recently emerged as the preferable method for volumetric imaging of the zebrafish brain. Based on z-scanning light-sheet illumination, the brain region can be excited slice-by-slice, and these side-sliced fluorescent images are built into the volumetric stack. This form of microscopy, originally designed for highly efficient optical sectioning, can effectively enhance the image contrast and reduce the phototoxicity in deep tissues [20].

According to Gaussian beam propagation [14], the relationship between the light-sheet thickness (2w0) and light-sheet nominal width (2z0) is given as where λ is the wavelength in vacuum, n is the refractive index of the medium on the object side, w0 is the waist radius of the Gaussian beam, and z0 is the Rayleigh length of the Gaussian beam. As illustrated, imaging with a wider FOV uses a thicker Gaussian beam.

Generally, thin light-sheet z-stacking images have desired optical sectioning and z-axis spatial resolution and do not require 3D deconvolution, since the thickness of light-sheet illumination is very close to the system depth of field (DOF) [20]. However, for zebrafish brain-wide imaging, the light sheet is thicker, where the fluorescence distributions outside the DOF will also be collected, blurring the image. Therefore, we resort to 3D deconvolution to redistribute the measured intensity from the z-stack images. Meanwhile, the nonuniform illumination of a Gaussian beam will produce an additional difference in each light-sheet image. In addition, melanogenesis tissues located at the surface of the fish head locally and significantly absorb the excitation energy, which also leads to dark stripes in the light-sheet image [24].

To address the abovementioned problems in zebrafish larva brain-wide imaging, we implement a thick Gaussian beam light-sheet microscope and expand the width of the light sheet to greater than 300 μm, which makes the thickness of the light sheet more than 8 μm, nearly 9.6-fold the system DOF. In this paper, we first derive an approximated 3D convolution forward model for thick light-sheet z-stacking imaging and propose a 3D deblurring method by introducing a Hessian regularization term to maintain the continuity of the neuron distribution. Employing this 3D reconstruction algorithm and a modified stripe-removal algorithm, the reconstructed z-stack images can exhibit high contrast and a high signal-to-noise ratio (SNR).

This paper is organized as follows. Section 2 introduces the forward model and 3D deblurring algorithm. Section 3 presents our experimental setup and preparation. Sections 4 and 5 present the results of both a numerical simulation and imaging experiments on fluorescent beads and a zebrafish larvae brain. Section 6 concludes and discusses the paper.

Assume that the illumination distribution ILS(x,y,z) is uniform in the x−y plane and symmetric with respect to the z-axis, i.e., ILS(x,y,z)=ILS(z)=ILS(−z). Assume that the system 3D point spread function (PSF) h(x,y,z) is space-invariant in the brain area. The imaging forward model of thick Gaussian beam light-sheet microscopy can be approximated as follows (see Appendix A): where f(x,y,z) is the true volumetric image, g(x,y,z=zj) is the 2D blurred image when the illumination and detection are located at the z=zj plane, hLS(x,y,z) is the overall 3D PSF affected by the illumination, and w=(p,q,r).

Here, we define a loss function as the sum of the fidelity term and Hessian regularization term, which is also popular in other imaging modalities [25,26]. The Hessian regularization term is from a priori knowledge of zebrafish brain (see Section 6). The fidelity term constrains the imaging forward model by using the mean square error, and the Hessian regularization term constrains the continuity of the neuron distribution. Such a loss function can be written as where α is the penalty parameter of the fidelity term. The Hessian regularization is defined as where αh and αz are the penalty parameters of continuity along the x−y plane and z-axis, respectively. The second-order partial derivatives of f in different directions are abbreviated as fi,i=xx,yy,zz,xy,xz,yz.

In this paper, we use the alternating direction method of multipliers (ADMM) [27] to solve this loss function based on its fast convergence performance in an L1-regularized optimization problem. The key of the algorithm is to decouple the L1 and L2 portions from the loss function [28]. Hence, we rewrite Eq. (3) by introducing auxiliary variables d and using the augmented Lagrangian methods (see Appendix B): where di,bi(i=xx,yy,zz,xy,xz,yz) are the auxiliary and dual variables, δ is the step size, and k is the iteration counter.

Consequently, the framework of the ADMM algorithm is presented in Algorithm 1.Table 1.

Evaluation Indicators of Different 3D Images in Fig. 5

In real experiments, light-sheet images are corrupted with photon shot noise and camera readout noise. We use the side window filtering (SWF) technique [29] to eliminate the mixed Poisson–Gaussian noise in the fluorescence image. In contrast to conventional local window-based filtering methods, the SWF method aligns the edges and/or corners of the window with the pixels being processed, which better preserves the image boundaries during the denoising process.

Furthermore, nonneuron melanogenesis tissues located at the surface of the fish head can absorb the excitation energy dramatically, which results in dark stripes along the excitation path. The widely used wavelet-FFT algorithm can be used to remove the stripe artifacts [30] but at cost price of deteriorating the image in regions devoid of stripes. Therefore, we introduce a prelocation method based on the output of the wavelet-FFT algorithm. Specifically, the positions and widths of all stripes are determined from the difference map between the original image and the wavelet-FFT processed image. Only the stripe areas in the original image will be multiplied with an adaptive Gaussian stripe function to correct for the artifacts. This modified algorithm can avoid information loss outside the stripe areas and provide consistent contrast across different regions of the image.

All larval zebrafish (elavl3:H2B-GCaMP6s, elavl3:GCaMP6s, and elavl3:EGFP) were raised in E3 embryo medium according to the standard protocol under 28.5°C. 0.2 mmol/L 1-phenyl-2-thiourea was added to the embryo medium to inhibit melanogenesis and allow optical imaging of the brain region. In our light-sheet imaging experiment, a 5–7 day post fertilization (dpf) zebrafish will be anesthetized before being embedded in a 1.5% low-melting agarose cylinder with a diameter of 1 mm and immersed in E3 medium. The brain area was electrically aligned to the center of the detection FOV being in good posture for z-stacking imaging.

The 3D PSF h(x,y,z) could be calibrated with sparse fluorescent beads statically embedded in the agarose cylinder or calculated by using the ImageJ plugin PSF generator and using the system parameters. The light-sheet illumination distribution ILS(z) could only be captured from the system utilizing imaging with a high-density fluorescent bead mixed medium. As shown in Eq. (2), the overall 3D PSF hLS(x,y,z) is given by the product of h(x,y,z) and ILS(z). After the imaging experiments, the results showed that the calculated PSF is better than the calibrated PSF during the reconstruction iteration (see Section 6).

The experiments are implemented on an HP Z6 Workstation (Intel Xeon Gold 6128 CPU @ 3.40 GHz × 24), the graphics card model is GeForce RTX 2080 Ti, and the operating system is Ubuntu 19.04. The 3D deconvolution algorithm is implemented using MATLAB R2019b.

Then, the ground truth was convolved with the overall 3D PSF function according to Eq. (2) followed by the introduction of Gaussian and Poisson mixed noise to generate a blurred 3D image stack [Fig. 5(b)]. These images were deconvolved with four deconvolution methods, the 2D Richardson–Lucy (RL) method, the 3D RL method [31,32], the 3D Wiener method [33], and our 3D method (also called the Hessian method) (see Section 2) [Figs. 5(c)–5(f)]. RL and Wiener are popular deconvolution algorithms for fluorescence microscopy images [34]. Of course, many software packages, such as DeconvolutionLab2 [35], also include the 3D RL and 3D Wiener algorithms and so on, but considering the consistency of the computing platform and the stability of performance, we choose the code provided by the MATLAB toolbox. We calculated four evaluation indicators, including the peak signal-to-noise ratio (PSNR) [36], SNR [33], structural similarity (SSIM) index [37], and correlation coefficient (R) [38] between the blurred image and four deblurred 3D image stacks, as shown in Table 1. The higher the indicator value, the closer the reconstructed image is to the real situation. Therefore, from this calculation, it is obvious that our method can recover more correct high-frequency components with much less noise.

There are three typical zebrafish lines frequently used for in vivo brain imaging, Tg (elavl3:EGFP), Tg (elavl3:GCaMP6s), and Tg (elavl3:H2B-GCaMP6s), which have neuronal-specific expression of GFP, the calcium indicator GCaMP6s localized in the cytoplasm of neurons, and the calcium indicator GCaMP6s localized in the nuclei of neurons, respectively.

For zebrafish larva brain-wide imaging, light-sheet illumination is required to cover a large FOV and possess concentrated energy in the focal plane. For this purpose, we implemented a thick Gaussian beam light-sheet microscope and expanded the width of the light sheet to more than 300 μm, making the thickness of the light sheet more than 8 μm, nearly 9.6-fold the system DOF. Our 3D deblurring method has been proposed to redistribute the measured intensity of each pixel in the light-sheet image to in situ voxels by 3D deconvolution. By introducing a Hessian regularization term to maintain the continuity of the neuron distribution and using a modified stripe-removal algorithm, the reconstructed z-stack images exhibit high contrast and a high SNR. These performance characteristics can facilitate subsequent processing, such as 3D neuron segmentation and recognition.

Comparing Fig. 5(c) with Fig. 5(f) in the simulation section, it has been verified that 3D deconvolution is more effective than 2D deconvolution for thick light-sheet imaging. In addition, comparing the different 3D deconvolution methods, especially referring to the reconstructed movies, our 3D deblurring method (see Fig. 3), including stripe-removal postprocessing, has been validated.

Here, we present some discussions concerning the deblurring methods and imaging settings. System.Xml.XmlElementSystem.Xml.XmlElementSystem.Xml.XmlElementSystem.Xml.XmlElementSystem.Xml.XmlElementSystem.Xml.XmlElementSystem.Xml.XmlElementSystem.Xml.XmlElement