The rat brain was prepared as reported previously (Pang et al., 2022). Briefly, adult animals were sacrificed by transcardial perfusion of 40 mL 1× phosphate buffered saline (PBS) and 40 mL 4% paraformaldehyde (PFA) successively. Subsequently, the sample preparation is followed by the Volumetric Imaging with Synchronized on-the-fly-scan and Readout (VISoR) imaging procedure (Wang et al., 2019). The sample was transferred into ice-cold 4% hydrogel monomer solution (HMS) for post-fixation at 4 degrees for 48 h. After post-fixation, the 4% HMS was replaced with a mixture of 15 mL 4% HMS and 15 mL 20% bovine serum albumin (BSA). The solution was degassed in a vacuum pump for 20 min with ice surrounding the centrifuge tubes. The sample was polymerized at 37°C for 4 h and rinsed with PBS three times to remove residual reagents. Next, the sample was sectioned into 300 μm thick slices. All slices were cleared with 4% sodium dodecyl sulfate (SDS) solution at 37°C for 24 h with gentle shaking. The slices were washed with PBS three times, 1 h for each round. The slices were immune-stained with anti-Aβ primary antibody (Biolegend No. 803002, 1:500 dilution with PBS) for 24 h at room temperature and washed with PBS for three times, 1 h for each round. The secondary antibody (JacksonImmuno Research No. 715-545-150, 1:200 dilution with PBS) was applied for 6 h at room temperature with gentle shaking and washed out with PBS for three times, 1 h per round. Finally, all slices were mounted on a customized slide (100 mm × 100 mm) for imaging (Figure 1A).

A modified VISoR microscope described in a previous study was used for all data collection (Xu et al., 2021). Two-channel imaging was carried out successively with 488 and 552 nm excitation. The emission light was collected through an Olympus 10 × 0.3 NA water immersion objective and filtered with bandpass filters (520/40, 600/50, from Semrock). Images were collected with a sCMOS camera (Flash 4.0 v3, Hamamatsu) (Figure 1A). All data were collected with a pixel resolution of 1 μm × 1 μm × 3.5 μm. The collected images were reconstructed into a dual channel whole rat brain for further segmentation and quantification (Figures 1B–D). The 3D reconstruction of the whole brain imaging dataset was the same as our previous study (Wang et al., 2019). Briefly, it mainly contains four steps. (i) Flattening the upper and lower surfaces of each brain section. (ii) Detecting the edges of adjacent sections and extracting the correspondences between two opposing surfaces. (iii) Morphing the correspondences of each section for limiting the morphological errors. (iv) Warping each section with the extracted and morphed correspondences.

A preprocessing pipeline of the dataset is presented to calibrate the brightness of serial brain sections as well as to enhance the signal-to-noise ratio. This step is crucial for the overall success of network training and testing. Two preprocessing strategies were applied successively:

Brightness calibration of the brain slices. The brightness of the imaging channel varies at different depths of the slices. Two reasons contribute most: (a) the difference in antibody concentration at different depths and (b) excitation light absorption by tissue through the optical propagation path. The brightness over z-axis was measured and automatically calibrated (Supplementary Figures 1A, B). Specifically, the corrected brightness of image stacks can be formulated by:

I(⋅) represent the intensity of each pixel. I_mean(⋅) indicates the mean intensity of each image slice. z means the slice index of the image stacks. The non-brain-slice pixels are excluded while calculating the mean intensity.

Signal-to-noise ratio (SNR) enhancement. As the Aβ plaques were labeled by immunofluorescent staining, the intensity of the fluorescence signal was relatively weak for image segmentation. The low SNR images will lead to poor performance in network training for image segmentation. To further enhance the segmentation performance, 3D Gaussian blur (σ_x,y,z = 2 for saving more details) was used to improve the SNR of all brain slices.

Due to the large size of 3D image stacks and the limitation of computer memory, we partitioned the image stacks into image blocks with size of 256 × 256 × 75. Image blocks having less than 1% brain tissue pixels were excluded for network training and testing.

Segmentation is a position-sensitive computer vision task. Due to the high density and small size of Aβ plaques in the whole rat brain, over-down-sampling of training data in the neural network will lead to the missing and position offset of small objects. The parallel structure of HRNet can combine high-resolution features with high-level semantic information. The high-resolution representations learned from the HRNet are not only spatially precise but also semantically strong.

The HRNet is connected in parallel (Figure 2A), which consists of parallel branches with high-to-low resolution. The resolution of the rth branch is 1/(2^r–1) of the resolution of the first stream, while the channel number is 2^r–1 of the first stream.

The basic module used in each subnetwork contains two 3 × 3 convolutional kernels where each kernel is followed by batch normalization (BN) and rectified linear unit (ReLU). Skip connection is used to connect the input and output of the module (Supplementary Figure 2A). The connection between different stages consists of transition modules and fusion modules (Supplementary Figure 2C), which are applied to exchange information between multi-resolution layers. Two output modes are provided (Supplementary Figure 2B): (i) HR-keep mode only outputs the high-resolution representation computed from the high-resolution convolution stream, (ii) HR-fuse mode combines the representations from all the high-to-low resolution parallel streams.

The training of neural networks is usually fully supervised and requires a large number of manual annotations, especially for image segmentation which requires pixel-level annotation. For 3D microscopic biomedical images, annotation requires extensive professional knowledge guidance, which undoubtedly leads to high expenses of money and time. In this paper, object-level weak annotations are used for pixel-level image segmentation tasks (Figure 2A). Specifically, we train the object detection framework with bounding box annotations and then segmentation is carried out by post-processing methods (Figure 2B and Supplementary Video 1).

Faster-RCNN (Ren et al., 2015) is a widely used object detection framework that usually consists of conv-body, region-proposal network (RPN), RoI-pooling, and classifier. In this work, a modified 3D form Faster-RCNN is used for detecting Aβ plaques in the volumetric image dataset. High-to-low resolution features of the 3D images are extracted by HRNet. After that, the features are input into the RPN and classifier to obtain the predicted bounding-box map. In Faster-RCNN, RPN depends on the default anchor settings, and we set the anchor size to fit the size range of Aβ plaques.

Subsequently, peak response mapping (PRM) (Zhou et al., 2018) is used to obtain additional visual cues to improve the segmentation performance (Figure 2A). The main idea of PRM is to generate a peak-response map by stimulating peaks in the class-aware map. Then, the most informative regions of each plaque are identified and mapped by the back-propagated peaks. Score maps generated from RPN in the Faster-RCNN framework can be regarded as a class of peak correspondence maps related to location (Dong et al., 2019). Therefore, we assume that peaks in the score map represent strong visual cues for the objects. The peak is back-propagated while the score map is sent to the classifier for further classification. The peak back-propagation can be interpreted as a random walk process. Each location in the bottom layer’s top-down relevance is formulated as its probability of being visited by the walker.

Consider a convolution layer with a filter size s×h×w. I_ijk and O_pqt are the spatial locations of the input and output feature maps, respectively. The visiting probability P(I_ijk) can be formulated by:

where the transition probability is formulated as:

I^i⁢j⁢k is the bottom-up activation value during the forward process. ReLU(W_{(i−p)(j−q)(k−t)}) means excluding negative weights as they are not helpful in improving the output response. Z_pqt is a factor to ensure ∑_p,q,t I_ijk|O_pqt = 1. Note that PRM only needs to be activated during testing rather than training.

Finally, we adopt an advanced 2D-OTSU (Zhang and Hu, 2008) algorithm for segmentation (Figure 2B). Unlike the original algorithm that uses handcrafted features as the second-dimension input, alternatively, we use the peak response map extracted from the neural network, which provides additional information in addition to grayscale intensity.

To quantitatively analyze the distribution of plaques in different brain regions, the 3D reconstructed rat brain needs to be registered to a rat brain atlas reference (Figure 2C). The WHS-SD rat brain atlas (Papp et al., 2014) is one of the most commonly used digital rat brain atlases. In this paper, we utilized atlas version 2 from public resources¹, which has 80 brain areas. We modified the images to remove the skull structures in the T2* template images, leaving only the brain structures.

The reconstructed 3D rat brain image dataset was registered to the WHS-SD brain atlas. Specifically, we used the WHS-SD rat brain template, annotation file, and atlas file provided by the public resources (see text footnote 1). The images were then converted to SimpleITK (Lowekamp et al., 2013) format to achieve fast alignment with the Elastix (Klein et al., 2009) toolbox. We adopted the non-rigid B-spline as the transform model. Mutual information (MI) and rigidity penalty was used as metrics for the similarity measure. The deformation field was optimized globally by using the stochastic gradient descent (SGD) algorithm. Order of B-spline interpolation was set to 3. Subsequently, a neuroanatomy expert was engaged to fine-tune the borderlines of all brain areas.

To comprehensively evaluate the accuracy of our weakly supervised segmentation method, we use three metrics for quantitative analysis, including the Dice score (DSC), Sensitivity (SST), and Hausdorff distance (HD). These metrics are defined as follows:

where TP denotes the true positive of the predicted pixels. FP denotes the false positive. FN denotes the false negative.

where X,Y denote the ground-truth pixel set and segmentation pixel set in the image segmentation task. sup represents the supremum. inf represents the infimum. d(⋅) is the distance metric.

Dice score is sensitive to interior pixels. SST represents the omission rate. HD is sensitive to boundary pixels. Note that 95% Hausdorff distance (HD95) is used here to remove the effect of minimal outliers.

We randomly selected 30 cropped blocks from the cerebral cortex as training data, and each block size was 256 × 256 × 75. Two experts were engaged to annotate images with bounding boxes. In addition, 28 blocks were randomly selected from different brain areas in three categories: cerebral cortex (Cortex), hippocampus (Hippo), and other brain areas (Other) as test data. All test image blocks were labeled with pixel-level annotations to evaluate the performance of our weakly supervised segmentation method. The preprocessing method was applied to each image block. All experiments were trained and tested with the PyTorch framework on a workstation with one NVIDIA Tesla V100S and 768 GB RAM. We carefully set the anchor sizes in range of [6,28] and the stride of RPN as 4. We used the SGD optimizer with the learning rate 0.01, the weight decay of 0.0001. The training process was terminated within 3,000 iterations.

We compared the proposed framework with widely used segmentation methods, including (a) traditional methods based on artificial features, such as Ilastik (Berg et al., 2019) and Segmenter-PMP34 (Chen et al., 2020); (b) fully supervised learning methods, such as U-Net (Ronneberger et al., 2015) and HRNet; and (c) weakly supervised learning methods with different weak labels, such as U-Net3D_rect and U-Net3D_grabcut (Khoreva et al., 2017). Different segmentation methods showed discernible results (Figures 3A–I). The evaluation metrics of different methods were calculated in distinct brain areas (Figures 3J–L and Supplementary Table 1).

Ilastik is a machine learning method based on the random forest algorithm. The user explicitly marks the features manually and applies batch processing to segment other images. However, this grayscale feature-based method showed the worst robustness through the whole brain, resulting in the poorest performance in the experiments. In addition, the computational cost of Ilastik is the highest due to CPU-based implementation. Segmenter_PMP34 is a conventional segmentation pipeline consists of Min-max intensity normalization, 2D Gaussian smoothing, 2D spot filter, watershed algorithm and size filter, which requires manual parameter adjustment while all methods in this experiment were fully automatic. For fairness, the parameters were adjusted in a cerebral cortex slice to reach the best performance and then applied to other brain slices. The metrics of the Hippo and other regions were significantly lower than our method. Interestingly, in the cerebral cortex region, the HD95 of Segmenter_PMP34 was better than our method but severely sacrificed SST indicating that many object pixels had not been detected.

To train the fully supervised network, two experts relabeled the training dataset with pixel-level labels. The U-Net architecture consists of an encoding and a decoding part. Encoding part repeatedly applies two 3 ×3 convolutional layers, each followed by a batch normalization layer and a ReLU. At each down-sampling step, the number of features is doubled. Decoding part recovers the original size by up-sampling the feature map. Every step of up-sampling consists of an up-sample layer that halves the number of features, and two 3 ×3 convolutional layers, each followed by a batch normalization layer and a ReLU. The final layer is a softmax layer. The HRNet architecture is the same as which used in our methods and HR-fuse mode is utilized for the best performance. We used the SGD optimizer with the learning rate 0.01, the weight decay of 0.0001 and momentum of 0.9. However, fully supervised segmentation requires a large number of pixel-level labeled image datasets. U-Net and HRNet produced poor segmentation results, partly because the training set was relatively small. All metrics were lower than our method, especially in SST.

We generated pixel-level labels through weakly labeled bounding boxes and then trained the semantic segmentation network iteratively (Khoreva et al., 2017). U-Net_3D_rect treats the pixels inside the 3D bounding box as pseudo labels. U-Net_3D_grabcut utilizes the GrabCut (Rother et al., 2004) algorithm inside the 3D bounding box to obtain the initial pseudo labels. We used these two pseudo labels to train the U-Net, respectively. The network architecture and parameter settings were the same as fully supervised network. The SST of these two methods was slightly better than our method, but Dice and HD95 were worse than ours as these methods treated too many background pixels as object pixels.

Furthermore, we clustered the plaques into three different categories according to their volumes, (a) small (<300 voxels), (b) medium (300–1,500 voxels), and (c) big (>1,500 voxels). The voxel size used here is 4 μm × 4 μm × 4 μm. We evaluated the performance of our method to analyze plaques of different sizes (Supplementary Figure 3 and Supplementary Tables 3, 4).

Our method only needed the bounding boxes as the weak labels. Almost all metrics achieved state-of-the-art performance. In addition, the ablation study was completed to verify the effectiveness of each module in our framework (Supplementary Table 2), which proved that HRNet (including fuse output mode), PRM, and 2D-OTSU all contributed to improving the segmentation performance.

To quantitatively analyze the distribution of Aβ plaques in the whole rat brain, segmented image blocks were montaged into brain slices and then reconstructed into the whole brain (Supplementary Videos 2, 3). Then we utilized a registration process to align the whole brain to the WHS-SD rat brain atlas (Figure 2C). The implementation of the registration method was based on SimpleITK (Lowekamp et al., 2013) and Elastix toolbox (Klein et al., 2009), performed on a workstation with 384 GB of RAM. The initial deformation field was subsequently obtained and fine-tuned by an expert. After registration, we applied the deformation field to the Aβ whole-brain segmentation binary masks. 3D rendering was performed for visualization of Aβ plaque distribution in the whole rat brain of several brain areas, including the cortex, hippocampus, and thalamus (Figures 4A–D). Aβ plaque maximum intensity projection in the coronal plane from the olfactory bulb to the caudal end revealed plaque density diversity in different brain areas or nuclei (Figure 4E). Finally, the Aβ plaque distribution was accurately calculated based on segmentation and registration by 61 brain regions that excluded several brain areas from the original atlas which did not contribute to our quantitative analysis, such as nerves, decussations, and commissures (Figure 5).

Most plaques were found in the neocortex, followed by the thalamus, brainstem, striatum, and other brain regions (Figure 5A). Due to the lack of a high-resolution rat brain atlas, we could not analyze our segmentation results by more fine-grained hierarchical brain areas or subregions. The total volume of segmented plaques in each region was similar to the tendency of total plaque counts (Figure 5B). We also measured the volume ratio of total plaque volume to the brain region volume by each brain region, and found that the volume ratio is higher in the frontal association cortex than in others (Figure 5C). We further calculated the plaque count density and the average plaque volume in each brain region (Figures 5D, E). Very few vascular-associated amyloid depositions were found in the early stage AD rat brain imaging data. Since we do not intend to distinguish these vascular-associated amyloid depositions from non-vascular amyloid depositions by our deep learning-based analysis pipeline, the segmentation result involves all kinds of amyloids.