Quantitative method of hypothalamic immunofluorescence image and system thereof

Abstract

This invention relates to the field of biomedical image processing technology, and discloses a method and system for quantitative analysis of hypothalamic immunofluorescence images. The method includes: performing spectral unmixing on multispectral fluorescence microscopy images based on a nonnegative matrix factorization algorithm to obtain a clean signal distribution map; using Gaussian Laplace filtering and watershed transform to achieve cell detection and segmentation; performing affine and B-spline registration between slice images and standard brain atlases to generate regions of interest masks for neural nuclei; using a local background adaptive correction strategy to perform fluorescence quantification and positive determination; and calculating Pearson correlation coefficient and Manders overlap coefficient to achieve colocalization analysis. The system includes a spectral unmixing module, a cell detection and segmentation module, an atlas registration and region recognition module, a fluorescence intensity quantification module, and a colocalization analysis and statistical output module.

Description

Technical Field

[0001] This invention relates to the field of biomedical image processing technology, and in particular to a quantitative method and system for hypothalamic immunofluorescence images based on multispectral fluorescence image analysis. Background Technology

[0002] Immunofluorescence imaging, an indispensable visualization tool in neuroscience and endocrinology research, uses specific fluorescently labeled antibodies to label and image target proteins in tissue sections. It has been widely applied to quantitative studies of neuronal distribution, protein expression levels, and intermolecular colocalization relationships in brain tissue. In studies of hypothalamic neuroendocrine regulation mechanisms, researchers typically need to simultaneously perform multi-channel immunofluorescence labeling of multiple biomarkers, such as gonadotropin-releasing hormone neurons, leptin receptors, and Kiss1 neuropeptide, within specific neural nuclei such as the arcuate nucleus and paraventricular nucleus, to reveal the spatial distribution characteristics and expression changes of these biomarkers under specific physiological and pathological conditions.

[0003] However, existing multichannel immunofluorescence image analysis methods face several technical bottlenecks in addressing the aforementioned research needs. Firstly, during multichannel fluorescence imaging, spectral crosstalk inevitably occurs between different fluorescent dyes due to partial overlap of their excitation and emission spectra. This results in interference components from other channel fluorescent dyes being mixed into the raw signals acquired from each channel. Traditional spectral separation methods primarily rely on linear regression or simple subtraction strategies, which are insufficient to fully eliminate multi-source crosstalk signals against complex tissue backgrounds. This is particularly problematic when the number of channels increases, leading to incomplete separation or signal distortion.

[0004] Chinese patent CN120927952B discloses a multiplex immunofluorescence detection method based on spatial reconstruction. Its core idea is to establish a three-dimensional voxel coordinate system using fluorescent reference beads for spectral unmixing and optimize it using three-dimensional spatial regularization constraints. This method achieves some success in reconstructing three-dimensional spatial continuity; however, its technical solution is mainly geared towards three-dimensional voxel reconstruction of continuous slices, focusing on the consistency of cross-slice spatial registration. It still falls short in terms of the efficiency and accuracy of refined spectral unmixing for two-dimensional single-slice multi-channel fluorescence images. Specifically, this method lacks a joint optimization mechanism for nonnegative matrix factorization and sparsity constraints of the spectral basis matrix, failing to achieve efficient adaptive crosstalk elimination within the two-dimensional image plane. Furthermore, this method does not involve automatic registration of tissue slice images with standard brain atlases to achieve automatic identification of regions of interest for specific neural nuclei. For small, poorly defined neural nuclei such as the arcuate nucleus and paraventricular nucleus of the hypothalamus, researchers still need to manually delineate the regions of interest, resulting in strong subjectivity and poor consistency.

[0005] Furthermore, in cell detection and counting, most laboratories still rely on manual visual counting or semi-automatic threshold segmentation for positive cell identification. This is not only time-consuming and labor-intensive, but also prone to missed or duplicate counts in densely populated cell areas, with significant differences in counting results between different operators. Regarding fluorescence intensity quantification, existing methods typically employ a uniform background subtraction strategy across the entire image, ignoring the spatial heterogeneity of autofluorescence intensity between different tissue regions. In brain tissue sections, due to the significant differences in autofluorescence levels between white and gray matter regions, a globally uniform subtraction strategy is prone to introducing systematic bias. In terms of colocalization analysis, conventional pixel-level overlap counting methods lack quantitative measures of spatial co-expression relationships of markers, making it difficult to accurately reflect, for example, the co-expression degree of leptin receptor and Kiss1 neuropeptide in the same neuronal population.

[0006] In summary, there is an urgent need for an intelligent quantitative method for immunofluorescence images that integrates functions such as spectral unmixing, automated cell detection, brain map registration, adaptive fluorescence quantification, and multi-channel colocalization analysis. This method would overcome problems such as incomplete elimination of spectral crosstalk, reliance on manual delineation of regions of interest, low cell counting efficiency, and inaccurate colocalization quantification in existing technologies. In particular, in animal model studies under special environmental conditions such as hypoxia and hypopressure, high-throughput and high-consistency quantitative analysis of large batches of slides is required, and the limitations of manual methods become increasingly apparent. Therefore, there is an urgent practical need to build an automated and standardized immunofluorescence image quantitative analysis pipeline. Summary of the Invention

[0007] To address the shortcomings of existing technologies, this invention provides a quantitative method and system for hypothalamic immunofluorescence images, aiming to solve the technical problems of difficulty in eliminating spectral crosstalk in multi-channel immunofluorescence images, difficulty in automatically identifying regions of interest in specific neural nuclei, low cell counting efficiency, and insufficient accuracy of colocalization analysis.

[0008] In a first aspect, this invention provides a method for quantitative analysis of hypothalamic immunofluorescence images, comprising the following steps: channel separation of acquired multispectral fluorescence microscopy images; spectral unmixing processing based on a nonnegative matrix factorization algorithm; applying nonnegative and sparsity constraints and iteratively updating to obtain a pure signal distribution map of each marker; automatic nucleus localization using Gaussian Laplace filtering and local maximum search on the cell nucleus staining channel images in the pure signal distribution map, combined with watershed transformation to achieve individualized segmentation of cells in dense regions; sequentially performing affine transformation and non-rigid B-spline deformation registration on the multispectral fluorescence microscopy images and standard brain atlases, mapping the boundaries of neural nuclei marked in the standard brain atlas to a pixel coordinate system to generate a region of interest mask; based on the region of interest mask and each cell mask region, performing modeling and subtraction of tissue autofluorescence using a local background adaptive estimation strategy, followed by fluorescence intensity quantification and positive cell counting; and performing colocalization analysis based on the fluorescence intensity of each marker channel after background correction and outputting the results of inter-group significant difference analysis.

[0009] Secondly, this invention provides a hypothalamic immunofluorescence image quantification system, including a spectral unmixing module, a cell detection and segmentation module, an atlas registration and region recognition module, a fluorescence intensity quantification module, and a colocalization analysis and statistical output module. Each module corresponds to a step in the aforementioned method. The spectral unmixing module is configured to perform spectral unmixing on multi-channel fluorescence images based on a non-negative matrix factorization algorithm to obtain a clean signal distribution map; the cell detection and segmentation module is configured to perform Gaussian Laplace filtering enhancement and watershed segmentation on the cell nucleus staining channel to output cell mask and centroid coordinates; the atlas registration and region recognition module is configured to perform affine and non-rigid registration between the slice image and a standard brain atlas to generate a region of interest mask for the neural nuclei; the fluorescence intensity quantification module is configured to perform fluorescence intensity calculation and positive cell determination using a local background adaptive correction strategy; and the colocalization analysis and statistical output module is configured to calculate the correlation and overlap coefficients between channels and perform inter-group statistical analysis.

[0010] The beneficial effects of this invention are as follows: High-precision spectral demixing of multi-channel fluorescence signals is achieved through a non-negative matrix factorization algorithm combined with sparsity constraints, effectively eliminating spectral crosstalk between fluorescent dyes; accurate identification of regions of interest in specific hypothalamic nuclei is achieved through automatic registration of standard brain maps, eliminating subjective biases from manual delineation; automatic individualized segmentation and counting of dense cells is achieved through a combination of Gaussian Laplace filtering and watershed transformation; spatial heterogeneity interference from tissue autofluorescence is eliminated through a local background adaptive correction strategy; and accurate quantification of spatial co-expression relationships of multiple markers is achieved through the joint calculation of Pearson correlation coefficient and Manders overlap coefficient. These steps form a deeply coupled data processing closed loop, where the region mask generated by map registration provides region-adaptive parameter feedback for spectral demixing, resulting in overall system performance significantly superior to the simple superposition of each step. Furthermore, the data flow between each step exhibits a cascading and closed-loop synergistic coupling relationship: the pure signal distribution map output from spectral demixing directly drives cell detection and segmentation; the spatial distribution information of cell centroid coordinates assists in feature point matching during spectral registration; and the region of interest mask generated by registration not only defines the spatial range for quantitative fluorescence but also feeds back the regional signal-to-noise ratio characteristics to the spectral demixing module to optimize its sparsity constraint parameters, thereby achieving a synergistic improvement in spectral demixing accuracy and regional adaptability. This deeply coupled multi-step architecture enables the present invention to achieve nonlinear gains that surpass the simple superposition of steps in terms of accuracy, automation, and repeatability in multi-channel immunofluorescence image quantitative analysis. Attached Figure Description

[0011] Figure 1 This is a schematic flowchart of the hypothalamic immunofluorescence image quantification method provided in an embodiment of the present invention.

[0012] Figure 2 This is a schematic diagram of the architecture of the hypothalamic immunofluorescence imaging quantitative system provided in an embodiment of the present invention. Detailed Implementation

[0013] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings. It should be understood that the following embodiments are only for explaining the present invention and do not constitute a limitation on the scope of protection of the present invention. Equivalent substitutions or improvements made by those skilled in the art based on the technical solutions disclosed in the present invention should all be covered within the scope of protection of the claims of the present invention.

[0014] See Figure 1 This invention provides a method for quantitative analysis of hypothalamic immunofluorescence images. The method comprises five tightly coupled processing steps, which form a closed-loop collaborative architecture through cascaded data flow and parameter feedback adjustment. Each step is described in detail below.

[0015] Step S1: Multichannel spectral demixing. This step aims to eliminate interchannel crosstalk caused by the overlap of fluorescent dye spectra during multichannel fluorescence imaging, and obtain the pure signal distribution map corresponding to each label, providing a reliable signal basis for subsequent cell detection and quantitative analysis.

[0016] In one embodiment of the present invention, the multispectral fluorescence micrograph to be processed is acquired by a fluorescence microscope equipped with a multi-band filter group. Preferably, the microscope objective has a magnification of 20x, a numerical aperture of 0.75, and an image pixel size of 0.325µm × 0.325µm. The number of fluorescence channels is set to 4 to 8 channels according to the experimental labeling scheme, and in a typical implementation scenario, 6 channels are set, corresponding to the DAPI nuclear staining channel, GnRH labeling channel, LepR labeling channel, Kiss1 labeling channel, and 2 auxiliary channels. The excitation wavelength and emission wavelength of each channel are configured according to the spectral characteristics of the selected fluorescent dye, wherein the excitation wavelength of the DAPI channel is 360nm and the emission wavelength is 460nm, the excitation wavelength of the Alexa Fluor 488 labeling channel is 488nm and the emission wavelength is 520nm, the excitation wavelength of the Cy3 labeling channel is 550nm and the emission wavelength is 570nm, and the excitation wavelength of the Alexa Fluor 647 labeling channel is 650nm and the emission wavelength is 670nm.

[0017] First, channel separation is performed on the original images of each channel, separating the multispectral image stack into independent single-channel two-dimensional images. Then, the mathematical model required for spectral demixing is constructed. Let the total number of channels be... The number of markers is The total number of pixels in the image is For the first image The pixel at the th position Observation signals on each channel This can be represented as a linear weighted superposition of the pure signals of each marker: ,in: For the first The first channel in the The observed gray value at each pixel ranges from 0 to 65535 (corresponding to a 16-bit image), and the unit is any fluorescence unit (au). The first in the spectral basis matrix Line 1 The element of the column represents the first element. The marker fluorescent dye in the first The relative emission intensity ratio on each detection channel, ranging from 0 to 1, is dimensionless. The first in the abundance matrix Line 1 The element of the column represents the first element. The marker in the first The pure signal strength at each pixel, taking values ​​in the range of non-negative real numbers, and the unit is au; This refers to the observation noise term, which includes photon shot noise and detector readout noise. (Subscript) Indicates the channel index. ; Subscript Indicates pixel index, ; Subscript Indicates the index of the marker. .

[0018] The above linear mixture model can be represented in matrix form as follows: ,in for The observation matrix for spectral basis matrix, for Abundance matrix. Spectral basis matrix. Reference samples stained with each fluorescent dye were collected before the formal experiment. Specifically, images of tissue sections labeled with only a single fluorescent dye were acquired in all channels, and the normalized average fluorescence intensity of each channel was filled into the corresponding column vector of the spectral basis matrix.

[0019] In obtaining the spectral basis matrix Subsequently, this invention employs a nonnegative matrix factorization algorithm to process the abundance matrix. Solve the problem. Preferably, based on the classical nonnegative matrix factorization objective function, introduce the L1 norm sparsity constraint to construct the following optimization objective function: ,in: is the Frobenius norm, representing the square root of the sum of squares of all elements in the matrix, used to measure the total residual between the observed matrix and the reconstructed matrix; is the L1 regularization coefficient, which physically controls the sparsity of the abundance matrix. The larger the value, the more likely the demixing result will compress weak signals to zero, thereby suppressing noise crosstalk. A smaller value retains more detailed signal but may leave behind crosstalk components; the value range is [value range missing]. to Preferably, the setting is adaptively configured based on the signal-to-noise ratio of each channel. Specifically, for channels with a signal-to-noise ratio below 10dB, Values For channels with a signal-to-noise ratio higher than 20dB, Values ; Abundance matrix elements The absolute value. Since a nonnegativity constraint has been applied. , here .

[0020] The above optimization problem is solved iteratively using a multiplicative update rule. In the... In this iteration, the update formula for the abundance matrix is: ,in: For the first The abundance matrix at the nth iteration Line 1 Column elements; For the first The updated value for the next iteration; For matrix The Line 1 The column elements represent the correlation intensity of the observed signal after projection onto the transpose of the spectral basis; For matrix The Line 1 The column elements represent the autocorrelation components of the current reconstructed signal; the transpose projection term in the numerator measures the strength of the observed data's association with each marker, the autocorrelation component in the denominator acts as a normalization, and the mathematical properties of the multiplicative update rule ensure that the initial value... During iteration, all elements remain non-negative, requiring no additional projection operations. The iteration terminates when the rate of change of the relative residual between two consecutive iterations is less than the convergence threshold. Or the number of iterations reaches the maximum value. Second-rate.

[0021] Preferably, the abundance matrix is ​​initialized using a non-negative random initialization strategy. Specifically, the abundance matrix is ​​initialized using a non-negative random initialization strategy. Each element is initialized to an interval Uniform random numbers within, where The abundance matrix is ​​the mean of all elements in the observation matrix. After the above iterative solution converges, the abundance matrix is... The The line is the first The pure signal intensity distribution of a marker across all pixels, rearranged into a two-dimensional image, yields the first... Pure signal distribution of the marker.

[0022] In one embodiment of the present invention, after the initial spectral demixing is completed in step S1, the region of interest mask from the spectral registration in step S3 can be fed back to this step. Specifically, the local signal-to-noise ratio statistics of pixels within each region of interest are used as region-level parameters, and the L1 regularization coefficient is adjusted accordingly. Perform regional adaptive adjustment—appropriately reduce the size of regions with strong fluorescence signals, such as the arcuate nucleus. To retain more signal details, the signal strength in areas with weaker signals, such as the periventricular nucleus edge, should be appropriately increased. To enhance noise suppression and thus achieve region-adaptive spectral unmixing optimization, this closed-loop feedback mechanism enables a balanced improvement in spectral unmixing accuracy across different tissue regions.

[0023] Step S2: Cell detection and individualized segmentation. This step takes the pure signal distribution map output from step S1 as input and performs automatic cell detection, localization and segmentation on the cell nuclear staining channel image to provide single-cell level spatial mask information for subsequent quantitative fluorescence and colocalization analysis.

[0024] First, the image corresponding to the DAPI nuclear staining channel is extracted from the clean signal distribution map obtained in step S1 as the source image for cell nucleus detection. In one embodiment of the present invention, the source image is first subjected to contrast normalization preprocessing, which linearly stretches the grayscale values ​​to a floating-point range of 0 to 1 to eliminate the influence of differences in imaging conditions between different slices on subsequent detection.

[0025] Subsequently, a Laplacian Gaussian filter was used to enhance the normalized nuclear staining image at multiple scales for speckle response. The core principle of the Laplacian Gaussian filter is to use parameters at different scales. The Laplacian response of the image is calculated below. Spherical or ellipsoidal nucleus structures will produce significant negative extremum responses at scales matching their size. Specifically, the formula for calculating the multi-scale Gaussian Laplacian response is: ,in: pixel coordinates In scale The filtered response value is a real number. The larger the absolute value of the negative response, the more likely there is a blob structure matching the current scale at that location. The standard deviation is A two-dimensional Gaussian function; and Gaussian functions along direction and Second-order partial derivatives in the direction; Represents a two-dimensional convolution operation; The input normalized kernel staining image at the pixel level The grayscale value at that location; This is the scale normalization factor, used to make response values ​​at different scales comparable. (Subscript) , These represent the horizontal and vertical pixel coordinates of the image, respectively.

[0026] Preferably, the scale parameter The target cell nucleus diameter range is set as a set of discrete values. In one embodiment of the present invention, the diameter of neuronal cell nuclei in hypothalamic tissue slices is typically between 5µm and 15µm, corresponding to approximately 15 to 46 pixels with a pixel size of 0.325µm. Therefore, Set as The filtering response is calculated at each of the five pixel scales, and the maximum value in the scale dimension is projected to obtain a comprehensive multi-scale response map.

[0027] After obtaining the multi-scale filtered response map, a local maximum search is performed to locate the coordinates of the candidate cell nuclei centers. Specifically, after taking a negative value for the response map (because the cell nucleus spots are negative extrema in the original response), a local maximum search is performed within a preset neighborhood radius. The search involves finding local maxima within the neighborhood. In one embodiment of the invention, the neighborhood radius... The value is set to half the average diameter of the target cell nucleus, approximately 5 pixels. A minimum response threshold for the maximum value is also set. To filter out false detection points caused by background noise, preferably, The standard deviation is set to twice the global standard deviation of the response map. After local maximum search, a set of candidate cell nucleus center coordinates is obtained. ,in The number of candidate cell nuclei detected.

[0028] In densely distributed regions, gray-level connectivity may exist between adjacent cell nuclei, necessitating watershed transformation for individualized segmentation. In one embodiment of this invention, before performing the watershed transformation, the filtered response map is first subjected to morphological H-minimum transformation processing, the parameters of which are... The dynamic range of the response map is preferably set to 5% to 10% to suppress oversegmentation caused by shallow minima. Then, the candidate nucleus center coordinates obtained in the previous step are applied as markers to the response map after H-minus transformation, and a marker-controlled watershed transformation is performed to segment the image into independent regions corresponding to each marker. Preferably, after segmentation, area constraint filtering is performed on each region to remove regions with areas smaller than the minimum nucleus area threshold (preferably 50 pixels, corresponding to approximately 5.3µm). 2 () or greater than the maximum cell nucleus area threshold (preferably 2000 pixels, corresponding to approximately 211µm) 2 The abnormal regions are identified, and the binary mask regions of each effective cell are ultimately output. and the corresponding centroid coordinates ,in This represents the number of effective cells after screening.

[0029] In one embodiment of the present invention, after completing the nuclear segmentation, the system further extends the segmentation results from the DAPI channel to other marker channels. Specifically, since the target protein of the immunofluorescence label may be located in different subcellular structures such as the nucleus, cytoplasm, or cell membrane, the present invention preferably extends the nuclear mask outward by a certain distance through a morphological dilation operation to generate an approximate mask for the cytoplasmic region. The dilation distance is preferably 3 to 8 pixels (corresponding to approximately 1.0 µm to 2.6 µm) to cover the fluorescence signal within the cytoplasmic range. For cell membrane markers, a larger dilation operation is preferably performed on the nuclear mask, with a dilation distance of 10 to 15 pixels (corresponding to approximately 3.3 µm to 4.9 µm), and the midline of the intersection with the dilated regions of adjacent cells is used as the basis for cell boundary division. The above-mentioned subcellular mask generation strategy ensures that the fluorescence intensity quantification of each marker channel corresponds spatially to the correct cell and subcellular region, avoiding the misattribution of signals from adjacent cells.

[0030] Step S3: Standard brain atlas registration and automatic region of interest (ROI) identification. This step is one of the core components with the highest technical barriers in this invention. It aims to automatically identify and delineate the boundaries of neural nuclei related to puberty development, such as the arcuate nucleus, paraventricular nucleus, and suprachiasmatic nucleus of the hypothalamus, by spatially registering tissue slice images with standard brain atlases, thereby generating a precise ROI mask. The output of this step not only provides spatial definition for quantitative fluorescence in step S4 but also provides regional parameter optimization information for spectral unmixing in step S1 through a feedback mechanism.

[0031] In one embodiment of this invention, the preset standard brain atlas is the Allen mouse brain reference atlas, which provides a set of whole-brain coronal cross-section reference images at a resolution of 25µm and fine boundary annotations of each neural nucleus. For the analysis of the hypothalamus region, this invention focuses on the arcuate nucleus (ARC), paraventricular nucleus (PVN), ventromedial nucleus (VMH), and suprachiasmatic nucleus (SCN), which are closely related to metabolic regulation and reproductive development.

[0032] The map registration process is divided into two stages, which are executed sequentially: global affine registration and local non-rigid fine registration.

[0033] In the global affine registration stage, the fluorescence microscopy images to be registered are first preprocessed. Specifically, the DAPI nuclear staining channel is extracted from the multi-channel image, or all channels are synthesized in grayscale to obtain a reference grayscale image containing tissue morphology information. Subsequently, the foreground mask of the tissue region is extracted using the Otsu automatic thresholding method, and morphological smoothing is performed on the mask boundaries to obtain tissue contour features. Simultaneously, the best-matching coronal section template image is selected from the standard brain atlas based on the anterior and posterior coordinates of the slice location as the registration target.

[0034] Preferably, the preliminary estimation of the slice location can be achieved through the following strategy: calculate the ratio of the tissue area of ​​the hypothalamus region to the cross-sectional area of ​​the third ventricle in the image to be registered, match this ratio with the theoretical ratio corresponding to each coronal section in the standard atlas, and select the section with the smallest matching error as the initial template. After determining the initial template, affine transformation registration is performed using it as a reference. The affine transformation includes 6 degrees of freedom parameters: horizontal translation... Vertical translation Rotation angle Horizontal scaling Vertical scaling and shear coefficient The objective function of the affine transformation is to maximize the normalized mutual information between the image to be registered and the template image. ,in: This is the normalized mutual information value, which takes a value greater than or equal to 1. The larger the value, the better the image quality. With template The stronger the statistical correlation between them; For the image to be registered The edge information entropy, measured in bits; template image Edge information entropy; The joint information entropy of the two images is given. This information entropy is obtained by performing joint histogram statistical calculation on the gray values ​​of the images, with the preferred gray-level bin number of the histogram being 64. The affine parameters are optimized using the Powell direction set optimization algorithm, with the goal of minimizing the negative value of the normalized mutual information. Preferably, a multi-resolution pyramid strategy is used to gradually refine the resolution from a coarse resolution of 4 times downsampling back to the original resolution to improve the convergence speed and avoid local extrema.

[0035] After global affine registration is completed, the local non-rigid fine registration stage begins. This stage employs a cubic B-spline deformation model to further correct local deformations of the tissue sections based on the affine registration results. Preferably, the spacing of the B-spline control point grid is adaptively set according to the image resolution and the degree of tissue deformation. In one embodiment, the initial control point spacing is set to 1 / 8 of the image's shorter side length, gradually refined to 1 / 16 in subsequent fine registration rounds. The objective function for non-rigid registration is: ,in, This is the overall objective function value for non-rigid registration. It takes the value of a real number, and the smaller the value, the higher the registration quality. For locally normalized cross-correlation coefficients, within a radius of... Calculation within the sliding window ( (Preferred to be 7 pixels), with a value range of -1 to 1. The larger the value, the higher the degree of local area matching. Indicates the image to be registered Through deformation field The transformed result image; Template image; The smoothness regularization term for the deformation field is defined as the sum of squares of the L2 norm of the gradient of the displacement vector at the control point. It is used to constrain the local smoothness of the deformation field to avoid unreasonable and drastic deformation. This is the regularization weight coefficient, and its value range is... to Preferred , The larger the value, the smoother the deformation, but the registration accuracy may decrease. Smaller mesh sizes allow for greater deformation but may lead to overfitting. The optimization process also employs a multi-resolution strategy, refining the mesh gradually in 3 to 4 levels from coarse to fine.

[0036] After registration, the labeled boundaries of each target neural nucleus in the standard brain atlas are mapped onto the pixel coordinate system of the original fluorescence microscopy image using an inverse deformation field, generating a region of interest mask for each nucleus. ,in The target number of nuclei. In one embodiment of the present invention, These correspond to the arcuate nucleus, paraventricular nucleus, ventromedial nucleus, and suprachiasmatic nucleus, respectively. It is worth noting that for small, poorly defined nuclei such as the arcuate nucleus and suprachiasmatic nucleus, the registered boundary may exhibit a positioning deviation of 1 to 3 pixels. Preferably, this invention performs a morphological dilation operation on the registered nucleus boundary to extend it by 2 pixels to ensure complete coverage of the target area. In subsequent quantitative fluorometry, background correction is used to eliminate non-specific signals introduced by the dilated region.

[0037] Furthermore, in one embodiment of the present invention, for small-volume nuclei such as the arcuate nucleus adjacent to the third ventricle, due to their significant morphological variation among different individuals and the lack of clear cellular architectural boundaries with the surrounding gray matter, the rigid boundary markings in the standard atlas may not accurately fit each experimental slice. Therefore, the present invention preferably introduces a boundary confidence assessment mechanism based on B-spline registration. Specifically, for each registered nucleus boundary, the deformation field gradient magnitude at the boundary pixel is calculated. When the gradient magnitude exceeds a preset threshold (preferably 15% of the control point spacing), the boundary segment is marked as a low-confidence region. In subsequent quantitative analysis, cells within the low-confidence region are assigned a weighting coefficient to reduce the impact of boundary uncertainty on the statistical results. In addition, the area, average gray value, and signal-to-noise ratio statistics of each region of interest obtained in this step are passed as feedback parameters to the spectral demixing module in step S1 for adjusting the L1 regularization coefficient. Regional adaptive adjustment is performed to form a closed-loop feedback path between step S3 and step S1.

[0038] Step S4: Fluorescence intensity quantification and positive cell determination. This step takes the pure signal distribution map output from step S1, the cell mask regions output from step S2, and the region of interest mask output from step S3 as inputs to perform background adaptive correction, fluorescence intensity quantification, and positive cell determination.

[0039] First, local background adaptive estimation is performed within each region of interest. Preferably, a sliding window percentile strategy is used to model the spatial distribution of tissue autofluorescence. Specifically, on the pure signal distribution map of each marker channel, a sliding window percentile strategy is employed, with a side length of... A square window of 1 pixel (preferred) (The number of pixels corresponds to a spatial range of approximately 16.25µm) is slid point by point within the region of interest, and the grayscale values ​​of the pixels within each window are sorted in ascending order. The first pixel is then selected. Percentile values ​​are used as local background estimates at the center of the window. Percentile parameters The preferred value is 10% to 25%, chosen because positive cells in immunofluorescence images typically occupy only a small portion of the region of interest. Therefore, a low percentile can better represent the autofluorescence level of nonspecific tissues without being artificially inflated by positive signals. After obtaining background estimates for all window center locations, two-dimensional cubic spline interpolation is used to fit the discrete background sampling points into a continuous background distribution surface. The background surface is subtracted pixel by pixel from the clean signal distribution map to obtain the background-corrected signal map.

[0040] In one embodiment of the present invention, considering that the autofluorescence levels of gray matter and white matter regions in brain tissue slices can differ by 2 to 5 times, a globally uniform background subtraction strategy may over-subtract weak positive signals in gray matter regions or under-subtract them in white matter regions. The aforementioned local adaptive strategy effectively addresses this heterogeneity through spatially varying background surfaces, ensuring that the quantitative baseline of positive signals in each region remains consistent.

[0041] After background correction, channel-level fluorescence intensity quantification was performed on each cell. For the cells detected in step S2... The number of effective cells, in the first Average fluorescence intensity on each marker channel With integrated optical density Calculated separately as follows:

[0042] ,

[0043] ,

[0044] in: For the first The cells in the first The average fluorescence intensity of each labeled channel, expressed in au; For the first Cell mask region The number of pixels contained; After spectral demixing in step S1, the first Each marker channel in pixel The pure signal strength at that location; For the first Each marker channel in pixel Local background fitting value at the location; The integral optical density value is equal to the sum of the background-corrected signals of all pixels within the mask area, reflecting the total expression level of the target marker in the cell.

[0045] Subsequently, a positive cell determination is performed. In one embodiment of the present invention, a positive determination threshold is set for each marker channel. The fluorescence intensity distribution was determined based on the negative control samples. Specifically, on negative control slides treated in the same batch as the experimental samples but without primary antibody, the fluorescence intensity distribution of all cells in the same region of interest was calculated at the [missing information - likely a specific point in the original text]. Average fluorescence intensity of each channel with standard deviation The positive threshold is set as follows: When the first Cell At that time, it was determined that the cell was in the first... Each marker channel was positive. The number of positive cells for each marker channel within each region of interest was counted. and positive rate .

[0046] Preferably, for cells located at the boundary of the region of interest, if the overlap area between the cell's mask region and the mask of the region of interest is less than 50% of the total cell area, the cell is excluded from the current nucleus count to avoid fluorescence intensity quantification deviation caused by boundary cutting. Further, in one embodiment of the present invention, normalized fluorescence intensity is also calculated for each cell. The original average fluorescence intensity is divided by the mean fluorescence intensity of the negative control in that channel to eliminate the influence of fluorescence efficiency differences between different channels on cross-channel comparisons. Furthermore, for multi-marker co-expression analysis, the system also generates multidimensional phenotypic vectors for each cell, encoding the positive or negative status of each channel as a binary vector. For example, for the GnRH, LepR, and Kiss1 channels, a cell's phenotypic vector of [1,0,1] indicates that the cell is GnRH positive, LepR negative, and Kiss1 positive. By statistically analyzing the phenotypic vectors of cells within each nucleus, the distribution of cell numbers with different phenotypic combinations can be obtained, providing a basis for the functional grouping of hypothalamic neurons.

[0047] Step S5: Colocalization analysis and statistical output. This step is based on the background-corrected fluorescence intensity data output in step S4. Colocalization analysis is performed on the channel pairs composed of the two selected marker channels. Combined with the quantitative results from the previous step, intergroup significance analysis is performed, and complete quantitative analysis results are output.

[0048] In the colocalization analysis phase, this invention employs two complementary colocalization indices: the Pearson correlation coefficient and the Manders overlap coefficient. For two specified marker channels... and In a region of interest Within, Pearson correlation coefficient The calculation formula is:

[0049] ,in: is the Pearson correlation coefficient, which ranges from -1 to 1. A positive value indicates that the two channel signals are positively correlated, i.e., co-location trend, while a negative value indicates a negative correlation, i.e., mutually exclusive distribution. Near 0 indicates no correlation. and Two channels at the pixel level Background-corrected fluorescence intensity at the location; and Two channels are located in the region of interest. The average fluorescence intensity within the range. The Pearson coefficient is sensitive to the overall linear correlation and can reflect the global covariance trend of the two markers in spatial distribution.

[0050] Manders overlap coefficients further quantify the degree of signal overlap between the two channels. For reference The formula for calculating the coefficient is:

[0051] ,in: For the channel The Manders overlap coefficient is used for reference, and its value ranges from 0 to 1. The larger the value, the more channels there are. A larger proportion of the signal is related to the channel. The positive areas overlap; This is an indicator function that takes the value 1 when the condition inside the parentheses is true and 0 otherwise. For channel The co-location threshold is preferably determined using the Costes automatic thresholding method. Specifically, this involves performing linear regression on a two-dimensional scatter plot of the two-channel signals and gradually decreasing the threshold until the regression intercept drops below zero; this point represents the threshold. Similarly, the threshold can be calculated using... For reference coefficient. and The asymmetry of coefficients can reveal the directionality of colocation relationships, for example... Gao Er Low indicates The markers are mainly distributed in within the positive area but Its distribution range is wider.

[0052] In one embodiment of the present invention, the colocalization heatmap is generated in the following manner: the normalized fluorescence intensity values ​​of the two channels are mapped to the horizontal and vertical axes of a two-dimensional scatter plot, respectively; two-dimensional kernel density estimation is performed on the scatter plot; and the density value is used as a color code to generate a heatmap. The high-density clustered areas on the heatmap intuitively reflect the colocalization intensity distribution pattern of the two markers.

[0053] In the statistical analysis output stage, this invention integrates the positive cell count, average fluorescence intensity, integrated optical density, Pearson coefficient, and Manders coefficient of each nucleus into a structured quantitative data table. For multi-group comparison scenarios, such as the comparison between the hypoxic-hypobaric group and the normoxic control group, the independent samples t-test or Mann-Whitney U test is used to assess the statistical significance of differences between groups, and Bonferroni correction or the Benjamini-Hochberg method is selected for multiple comparison correction according to the number of groups. Significance level The preferred value is set as follows .

[0054] Preferably, before performing statistical tests, the system first performs a normality test on the quantitative data of each group. Specifically, the Shapiro-Wilk test is used to determine the normality of the positive cell count and average fluorescence intensity data of each group. Normality is determined when the sample size is greater than or equal to 3 and... If the distribution is determined to be normal, a parametric test, namely the independent samples t-test, is used; otherwise, a non-parametric test, namely the Mann-Whitney U test, is used. For comparisons involving three or more groups, this invention uses one-way ANOVA or the Kruskal-Wallis test to assess overall differences, and performs post-hoc tests on pairwise comparisons only if the overall differences are significant.

[0055] In terms of report output, this invention generates multi-dimensional visual output content: grouped bar charts are generated for the positive cell counts of each neural nucleus, where error bars represent standard error or standard deviation; grouped box plots are generated for the average fluorescence intensity of each marker channel to show the data distribution characteristics; a two-dimensional scatter density heatmap is generated for the colocalization analysis results; and radar charts are generated for the changes in quantitative indicators of each nucleus under different experimental conditions to intuitively display the comprehensive difference patterns of multi-dimensional characteristics. The final output includes a complete quantitative analysis report containing all the above charts and statistical test results, providing precise and objective image-based quantitative evidence for elucidating the molecular mechanism by which obesity affects the initiation of puberty in high-altitude environments through the hypothalamic pathway.

[0056] See Figure 2 The present invention also provides a hypothalamic immunofluorescence image quantification system, which includes five functional modules, each corresponding to one of the five steps in the above method embodiments, together forming a complete image quantification analysis pipeline.

[0057] The spectral unmixing module, as the system's primary processing unit, is configured to receive multi-channel raw image data acquired by fluorescence microscopy and perform channel separation and spectral unmixing operations based on nonnegative matrix factorization. Internally, this module includes a spectral basis matrix storage subunit and an iterative solution engine subunit. The spectral basis matrix storage subunit maintains single-stain reference spectral data corresponding to each experimental scheme, while the iterative solution engine subunit performs iterative calculations of the multiplicative update rule and receives region-level signal-to-noise ratio parameters from the map registration and region recognition module to adaptively adjust the L1 regularization coefficient. This module outputs the pure signal distribution maps of each marker and transmits them to the downstream cell detection and segmentation module and fluorescence intensity quantification module.

[0058] The cell detection and segmentation module is configured to perform automatic cell detection and segmentation on the DAPI nuclear staining channel in the clean signal distribution map output by the spectral unmixing module. This module includes a multi-scale filtering subunit, a maximum detection subunit, and a watershed segmentation subunit. The multi-scale filtering subunit performs Gaussian Laplacian convolution operations on the nuclear staining image sequentially according to a preset scale parameter sequence and outputs a multi-scale maximum response map. The maximum detection subunit locates the coordinates of the candidate cell nuclei's centers on the response map through neighborhood search. The watershed segmentation subunit performs a marker-controlled watershed transformation using the candidate center coordinates as markers and filters out abnormal regions by area selection. The module ultimately outputs the binary mask region and centroid coordinates of each effective cell.

[0059] The map registration and region identification module is the most technically complex core module in this system, configured to spatially register fluorescence microscopy images with pre-stored standard brain maps. This module includes a template matching subunit, an affine registration subunit, a non-rigid registration subunit, and a nucleus annotation mapping subunit. The template matching subunit automatically selects the best-matching coronal section template from the standard map library based on tissue morphology characteristics; the affine registration subunit calculates global affine transformation parameters by optimizing normalized mutual information; the non-rigid registration subunit performs fine registration using a multi-resolution B-spline deformation model based on the affine results; and the nucleus annotation mapping subunit maps the nucleus boundaries in the map to the original image coordinate system using an inverse deformation field to generate a region of interest mask. Preferably, this module also passes the local signal-to-noise ratio characteristics of each region of interest as feedback parameters to the spectral unmixing module, forming a closed-loop collaborative optimization between unmixing accuracy and region adaptability.

[0060] The fluorescence intensity quantification module is configured to perform local background adaptive estimation, fluorescence intensity quantification, and positive cell determination under the dual constraints of region-of-interest (ROI) mask and cell mask. This module includes a background modeling subunit, an intensity calculation subunit, and a positive cell determination subunit. The background modeling subunit uses a sliding window percentile strategy to generate a continuous background distribution surface; the intensity calculation subunit calculates the average fluorescence intensity and integrated optical density of each channel on a cell-by-cell basis after background subtraction; the positive cell determination subunit determines a threshold based on statistical parameters of the negative control and performs binarization determination.

[0061] The colocalization analysis and statistical output module is configured to receive output data from the fluorescence intensity quantification module, calculate the Pearson correlation coefficient and Manders overlap coefficient of the channel pair composed of the two selected marker channels, generate a colocalization heatmap, and integrate quantitative indicators from each step to perform statistical tests on inter-group differences. This module outputs all analysis results in the form of structured data tables and multidimensional visualization charts, forming a complete quantitative analysis report.

[0062] In one embodiment of the present invention, the above five modules are deployed on a computing platform equipped with a multi-core processor and a graphics acceleration unit, and the modules exchange information through a standardized data interface. Preferably, the computing platform is configured with an Intel Xeon series or AMD EPYC series processor (no less than 8 cores), equipped with no less than 32GB of system memory to meet the memory residency requirements of large-size image data, and the graphics acceleration unit is preferably a GPU supporting CUDA parallel computing (with no less than 8GB of video memory) to accelerate matrix operations for spectral unmixing and interpolation calculations for B-spline deformation fields. The system software environment is built on Python 3.8 or above, and the core dependent libraries include NumPy and SciPy for matrix operations and optimization solutions, scikit-image for image filtering and segmentation, SimpleITK or ANTsPy for spectral registration, and matplotlib for visualization output.

[0063] When processing a multispectral fluorescence microscopy image containing 6 channels and a resolution of 2048×2048 pixels, the typical processing time distribution for each module is as follows: spectral unmixing module takes approximately 8 seconds (of which iterative solution takes approximately 6 seconds), cell detection and segmentation module takes approximately 5 seconds, map registration and region recognition module takes approximately 12 seconds (of which B-spline registration takes approximately 9 seconds), fluorescence intensity quantification module takes approximately 3 seconds, and colocalization analysis and statistical output module takes approximately 2 seconds, with a total processing time of approximately 30 seconds. Compared to researchers manually completing the nucleus delineation, cell counting, and fluorescence quantification of the same image, which typically takes 60 to 90 minutes, this system achieves a processing speed improvement of over 100 times. The system provides batch processing functionality, supporting automated pipeline processing of multiple slice images from the same experimental batch, and automatically summarizing the quantitative results of each slice into an integrated analysis dataset. Data exchange between modules uses a unified data structure encapsulation to ensure the integrity and traceability of information in the processing chain. Specifically, the data objects passed between modules comprise three layers: image data bodies, metadata tags, and processing logs. The image data bodies carry the actual image matrix and calculation results; the metadata tags record information such as image acquisition parameters, channel configurations, pixel dimensions, and coordinate systems to ensure consistent data interpretation across modules; and the processing logs automatically record parameter settings, iteration processes, and quality control indicators for each module, facilitating retrospective review and problem localization. The system also provides a parameter configuration file interface, allowing researchers to flexibly adjust key parameters of each module according to different experimental protocols and organizational types without modifying the core algorithm code, thus ensuring the system's universal adaptability and scalability across various experimental scenarios.

[0064] To verify the effectiveness of the method of this invention, a systematic performance test was conducted on the hypothalamic immunofluorescence images of a hypoxic-hypobaric obese mouse model. The test dataset contained 72 hypothalamic coronal section immunofluorescence images from 6 mice in the hypoxic-hypobaric obese group and 6 mice in the normoxic normal weight control group. Six sections were collected from each mouse at different anterior and posterior coordinate positions. Each section contained six fluorescence channels: DAPI, GnRH, LepR, Kiss1, NeuN, and GFAP. The image resolution was 2048×2048 pixels.

[0065] Regarding spectral unmixing performance, the evaluation was conducted using a standard validation sample with a known mixing ratio of manually prepared fluorescent dyes as a benchmark. The non-negative matrix factorization (NMF) spectral unmixing method of this invention achieved an average crosstalk residual rate of 2.3% in a 6-channel configuration, a reduction of approximately 73.6% compared to the 8.7% crosstalk residual rate of the traditional linear regression method, and a reduction of approximately 55.8% compared to the three-dimensional voxel space regularized unmixing method in the comparative document (which had a crosstalk residual rate of 5.2% in a two-dimensional single-slice scenario). The introduction of sparsity constraints resulted in an average improvement of 4.2 dB in the signal-to-noise ratio of each channel after unmixing compared to the standard NMF method without sparsity constraints.

[0066] Regarding cell detection accuracy, manual counting by experienced researchers on the same images was used as the gold standard. In the arcuate nucleus region, the cell detection precision of the method of this invention was 94.7%, the recall was 91.3%, and the F1 score was 93.0%. In the paraventricular nucleus region, where cell density is highest, the watershed transformation effectively solved the cell adhesion problem, and the F1 score remained above 89.5%, while the F1 score of the traditional global threshold segmentation method in this region was only 72.8%.

[0067] Regarding the accuracy of atlas registration, the nucleus boundaries manually annotated by professional neuroanatomists were used as the reference standard. The affine-plus-B-spline combined registration method of this invention achieved a Dice similarity coefficient of 0.87 in the arcuate nucleus region, 0.91 in the paraventricular nucleus region, and 0.83 in the suprachiasmatic nucleus region. Compared to the method using only affine registration, the non-rigid B-spline fine registration improved the Dice coefficient by an average of 0.09. The average localization error of the registered nucleus boundaries was controlled within 1.5 pixels (approximately 0.49µm), meeting the accuracy requirements for quantitative analysis at the hypothalamic neural nucleus level.

[0068] Regarding the accuracy of quantitative fluorescence, the local background adaptive correction strategy, compared to the global mean background subtraction method, reduced the coefficient of variation of the average fluorescence intensity of positive cells from 18.3% to 9.6%, effectively eliminating the systematic bias caused by the difference in autofluorescence between gray and white matter regions. The batch-to-batch coefficient of variation of the integrated optical density value decreased from 14.7% to 7.2%, indicating that this strategy significantly improved the batch-to-batch reproducibility of the quantitative results.

[0069] In terms of colocalization analysis, the Pearson correlation coefficient calculated for the GnRH and Kiss1 channel pairs was 0.72, and the Manders M1 coefficient was 0.81 (with GnRH as a reference). These quantitative indicators are consistent with the high co-expression relationship between Kiss1 neurons and GnRH neurons in the arcuate nucleus region reported in the literature. The difference in GnRH-positive cell counts in the arcuate nucleus region between the hypoxic-hypotensive obesity group and the control group reached a statistically significant level. Kiss1 positive cell counts also showed significant differences between groups ( The differences between groups of markers in the paraventricular nucleus region were not statistically significant. These results are consistent with the arcuate nucleus-specific hypothesis that obesity under hypoxic and hypobaric conditions affects the initiation of puberty, and verify that the method of this invention has sufficient sensitivity and reliability in detecting such region-specific biological effects.

[0070] In summary, the method of the present invention exhibits significantly superior performance compared to existing technologies in all dimensions, including spectral unmixing accuracy, cell detection accuracy, map registration accuracy, fluorescence quantitative repeatability, and colocalization analysis sensitivity. The deep coupling and closed-loop feedback mechanism of the five processing steps enable the overall system performance to achieve synergistic gains that surpass the simple summation of the individual steps.

[0071] It is particularly noteworthy that the contribution of the closed-loop mechanism of feeding back the region of interest mask generated by the spectral registration in step S3 to the spectral demixing module in step S1 to the overall system performance has been verified by a specific ablation experiment. With the region-adaptive regularization parameter feedback enabled, the residual rate of spectral demixing crosstalk in the arcuate nucleus region further decreased from 3.1% without feedback to 2.3%, a reduction of 25.8%. Simultaneously, in small nucleus regions with weak signals, such as the suprachiasmatic nucleus, the closed-loop feedback increased the recall rate of positive cell detection from 85.2% to 91.3%, confirming the signal preservation effect of region-adaptive sparsity constraints on weak signal regions. The above experimental results demonstrate that the deeply coupled closed-loop architecture between each step produces significant synergistic gains, with the overall performance improvement exceeding the upper limit achievable by independent optimization of each step, fully reflecting the nonlinear system advantage of this invention where one plus one equals more than two.

[0072] The embodiments of the present invention are not limited to the specific embodiments described above. Those skilled in the art can make various equivalent changes or substitutions based on the technical solutions of the present invention, and all such changes or substitutions should be included within the protection scope of the present invention.

Claims

1. A method for quantifying hypothalamic immunofluorescence images, characterized by, The method comprises the following steps: Channel separation is performed on the collected multi-spectral fluorescence microscopic image, spectral unmixing processing is performed based on a non-negative matrix factorization algorithm, an observation matrix is constructed based on each channel observation signal, a spectral basis matrix is constructed based on each fluorescent dye reference emission spectrum, a residual error of a product of the observation matrix, the spectral basis matrix and an abundance matrix is minimized and solved, and a non-negative constraint and a sparsity constraint are applied, and a pure signal distribution map of each marker is obtained through iterative updating; Cell detection and segmentation are performed on a cell nucleus staining channel image in the pure signal distribution map, multi-scale spot response enhancement is performed by using a Gaussian Laplace filter, local maximum value search is performed on a filter response map to locate candidate cell nucleus center coordinates, and dense region cell individual segmentation is realized by using the candidate cell nucleus center coordinates as seed points through a watershed transform, and each cell mask region and a centroid coordinate are output; Atlas registration is performed on the multi-spectral fluorescence microscopic image and a corresponding coronal section template image in a standard brain atlas, tissue contour features and gray distribution features are extracted, and affine transformation and non-rigid B-spline deformation registration are sequentially performed with reference to the corresponding coronal section template image in the standard brain atlas, a neural nucleus boundary labeled in the standard brain atlas is mapped to a pixel coordinate system, and a region of interest mask of each target neural nucleus is generated; Based on the region of interest mask and each cell mask region, a local background adaptive estimation strategy is used to model and subtract tissue autofluorescence, the average fluorescence intensity and integral optical density value of each cell in each marker channel after background correction are calculated, positive determination is performed according to a positive determination threshold, and a positive cell count of each nucleus is counted; Based on the fluorescence intensity of each marker channel after background correction, co-localization analysis is performed on a channel pair composed of two selected marker channels, a correlation coefficient and an overlap coefficient are calculated, a co-localization heat map is generated, significant difference analysis between groups is performed in combination with the positive cell count and the fluorescence intensity statistical value, and a quantitative analysis result is output.

2. The method of claim 1, wherein, The number of acquisition channels of the multi-spectral fluorescence microscopic image is 4 to 8, the cell nucleus staining channel adopts a channel with an excitation wavelength of 360 nm and an emission wavelength of 460 nm corresponding to a DAPI dye, and the iteration termination condition of the non-negative matrix factorization algorithm is that the residual error change rate is less than a preset convergence threshold or the number of iterations reaches a preset maximum number of iterations.

3. The method of claim 1, wherein, The scale parameter of the Gaussian Laplace filter is set to a plurality of discrete values corresponding to the target cell nucleus diameter range, the neighborhood radius of the local maximum value search is determined according to the average value of the target cell nucleus diameter, and a morphological H-minima transform is performed on the filter response map before the watershed transform to suppress over-segmentation.

4. The method of claim 1, wherein, The standard brain atlas is an Allen mouse brain atlas or a Paxinos mouse brain atlas, the parameters of the affine transformation include a translation vector, a rotation angle and an anisotropic scaling ratio, and the control point grid spacing of the non-rigid B-spline deformation registration is adaptively set according to the image resolution and the tissue deformation degree.

5. The method of claim 1, wherein, The sparsity constraint in the non-negative matrix factorization algorithm is achieved by introducing an L1 norm regularization term of the abundance matrix into the objective function. The regularization coefficient of the L1 norm regularization term is adaptively adjusted according to the signal-to-noise ratio of each channel, and the iterative update adopts a multiplicative update rule to ensure that each element of the abundance matrix always satisfies the non-negativity constraint during the iteration process.

6. The method of claim 1, wherein, In the spectral registration, the affine transformation obtains the optimal affine parameters by maximizing the normalized mutual information between the multispectral fluorescence microscopy image and the coronal cross-section template image. The non-rigid B-spline deformation registration uses the output of the affine transformation as the initial value and optimizes the displacement of the B-spline control points by minimizing the local normalized cross-correlation loss function. During the optimization process, a deformation field smoothness regularization term is introduced to constrain the deformation gradient.

7. The method of claim 1, wherein, The local background adaptive estimation strategy includes: traversing the image in each region of interest using a preset sliding window, sorting the pixel gray values ​​in each window, taking the lower percentile value as the local background estimate at the center of the window, generating a continuous background distribution surface by interpolating and fitting the background estimates at the centers of all windows, and subtracting the background distribution surface from the clean signal distribution map.

8. The method of claim 1, wherein, In the colocation analysis, the correlation coefficient is the Pearson correlation coefficient, the overlap coefficient is the Manders overlap coefficient, and the Manders overlap coefficient includes the M1 coefficient and the M2 coefficient, which are respectively referenced by the two target channels. The colocation heatmap is generated by mapping the normalized intensity values ​​of the two target channels to the horizontal and vertical axes of a two-dimensional scatter plot and then performing kernel density estimation.

9. The method of claim 1, wherein, The clean signal distribution map of the spectral unmixing process is constrained by the region of interest mask generated by the spectral registration, and the local signal-to-noise ratio characteristics of each region of interest are fed back to the sparsity constraint parameters of the non-negative matrix factorization algorithm to achieve regional adaptive spectral unmixing optimization.

10. A hypothalamic immunofluorescence image quantification system for implementing the method according to any one of claims 1 to 9, characterized in that, include: The spectral unmixing module is configured to perform channel separation on the acquired multispectral fluorescence microscopy images. Based on the non-negative matrix factorization algorithm, it performs spectral unmixing processing on the separated multi-channel images. By minimizing the residual between the product of the observation matrix, the spectral basis matrix, and the abundance matrix, and applying non-negative and sparsity constraints, it iteratively updates the pure signal distribution map corresponding to each marker. The cell detection and segmentation module is configured to perform Gaussian Laplacian filtering enhancement and local maximum search localization on the cell nucleus staining channel image in the pure signal distribution map, and perform watershed transformation based on the coordinates of the candidate cell nucleus center to achieve individualized cell segmentation, and output the mask region and centroid coordinates of each cell; The map registration and region recognition module is configured to sequentially perform affine transformation and non-rigid B-spline deformation registration on the multispectral fluorescence microscopy image and the preset standard brain map, map the boundaries of the neural nuclei marked in the standard brain map to the pixel coordinate system, and generate a region of interest mask corresponding to each target neural nucleus. The fluorescence intensity quantification module is configured to model and subtract tissue autofluorescence based on the region of interest mask and the mask region of each cell using a local background adaptive estimation strategy, calculate the average fluorescence intensity and integrated optical density value of each cell in each marker channel, and perform positive determination and positive cell count statistics. The colocalization analysis and statistical output module is configured to calculate the Pearson correlation coefficient and Manders overlap coefficient between channel pairs consisting of two selected marker channels based on the fluorescence intensity of each marker channel after background correction, generate a colocalization heatmap, and perform intergroup significance analysis by combining positive cell count and fluorescence intensity statistics, and output quantitative analysis results.

Images