Advertisement

Mouse Mammary Gland Whole Mount Density Assessment across Different Morphologies Using a Bifurcated Program for Image Processing

Published:September 13, 2022DOI:https://doi.org/10.1016/j.ajpath.2022.06.013
      Mammographic density is associated with increased breast cancer risk. Conventional visual assessment of murine mouse models does not include quantified total density analysis. A bifurcated method was sufficient to obtain relative density scores on a broad range of two-dimensional whole mount images that contained both normal and abnormal findings. Image processing techniques, including a ridge operator and a gaussian denoising method, were used to isolate background away from mammary epithelium and use mean pixel intensity to represent mammary density on genetically engineered mouse models for breast cancer in mice 4 to 29 months of age. The bifurcated method allowed for application of an optimal image processing approach for the structural elements present in the whole mount images. Gaussian denoising was the optimal approach when more dense lobular growth and tertiary branching dominate and a ridge operator when epithelial growth was more sparse and secondary branching was the more dominant structural feature. The two processing approaches were combined in a single experimental flow program using an initial image density measurement as the decision point between the two approaches. Higher density was associated with lobular growth, tertiary branching, fibrotic stroma, and presence of cancer. The significance of the study is development of a readily accessible program for digital assessment of mammary gland whole mount density across a range of mammary gland morphologies.

      Graphical abstract

      In the United States, the primary means of screening for breast cancer is the mammogram, which is used to evaluate both the presence of radio-dense lesions and overall breast density.
      • Fuller M.S.
      • Lee C.I.
      • Elmore J.G.
      Breast cancer screening: an evidence-based update.
      Effective use of population-based mammography screening is limited to women older than 40 years of age because at younger ages hormonal-related breast density can result in either false-positive results from what are actually normal premenopausal growth patterns or false-negative results due to pathologic lesions being indistinguishable from normal premenopausal breast.
      • Burton A.
      • Maskarinec G.
      • Perez-Gomez B.
      • et al.
      Mammographic density and ageing: a collaborative pooled analysis of cross-sectional data from 22 countries worldwide.
      ,
      • Carney P.A.
      • Miglioretti D.L.
      • Yankaskas B.C.
      • et al.
      Individual and combined effects of age, breast density, and hormone replacement therapy use on the accuracy of screening mammography.
      A second use of mammography is to assess breast cancer risk
      • Mendes J.
      • Matela N.
      Breast cancer risk assessment: a review on mammography-based approaches.
      because the percentage of mammographic density is positively associated with breast cancer risk.
      • Boyd N.F.
      • Martin L.J.
      • Yaffe M.J.
      • Minkin S.
      Mammographic density and breast cancer risk: current understanding and future prospects.
      • Jakes R.
      • Duffy S.
      • Ng F.
      • Gao F.
      • Ng E.
      Mammographic parenchymal patterns and risk of breast cancer at and after a prevalence screen in Singaporean women.
      • Huo Z.
      • Giger M.L.
      • Wolverton D.E.
      • Zhong W.
      • Cumming S.
      • Olopade O.I.
      Computerized analysis of mammographic parenchymal patterns for breast cancer risk assessment: feature selection.
      In murine models, mammary gland whole mounts serve as the standard evaluation of whole mammary gland for abnormal growth patterns and preneoplastic lesions. Conventional visual assessment includes evaluation of branching structures, hyperplastic alveolar nodules, and lobular growth.
      • Dabydeen S.A.
      • Kang K.
      • Díaz-Cruz E.S.
      • Alamri A.
      • Axelrod M.L.
      • Bouker K.B.
      • Al-Kharboosh R.
      • Clarke R.
      • Hennighausen L.
      • Furth P.A.
      Comparison of tamoxifen and letrozole response in mammary preneoplasia of ER and aromatase overexpressing mice defines an immune-associated gene signature linked to tamoxifen resistance.
      • Alothman S.J.
      • Wang W.
      • Goerlitz D.S.
      • Islam M.
      • Zhong X.
      • Kishore A.
      • Azhar R.I.
      • Kallakury B.V.
      • Furth P.A.
      Responsiveness of Brca1 and Trp53 deficiency–induced mammary preneoplasia to selective estrogen modulators versus an aromatase inhibitor in Mus musculus.
      • Nakles R.E.
      • Kallakury B.V.S.
      • Furth P.A.
      The PPARγ agonist efatutazone increases the spectrum of well-differentiated mammary cancer subtypes initiated by loss of full-length BRCA1 in association with TP53 haploinsufficiency.
      Approaches for quantitative analysis of mammary gland structure using digitized images have been described as time-consuming and sometimes frustrated by uneven or noisy images.
      • Fernandez-Gonzalez R.
      • Barcellos-Hoff M.H.
      • Ortiz-de-Solórzano C.
      Quantitative image analysis in mammary gland biology.
      To address this need, previously described methods have been published with effective results using Sholl analysis in ImageJ software version 1.49c (NIH, Bethesda, MD; http://imagej.nih.gov/ij) for branching density and a Z-Stack image processing protocol.
      • Stanko J.P.
      • Easterling M.R.
      • Fenton S.E.
      Application of Sholl analysis to quantify changes in growth and development in rat mammary gland whole mounts.
      ,
      • McGinley J.N.
      • Thompson H.J.
      Quantitative assessment of mammary gland density in rodents using digital image analysis.
      However, there remains a need for a protocol that can be efficiently applied to conventional two-dimensional images and represent total epithelial density. Moreover, the mammary gland is not a static structure in adult mice but rather undergoes cyclic changes with each estrous cycle, with lobuloalveolar budding dominating during diestrus and more simple secondary branching patterns during estrus.
      • Snijders A.M.
      • Langley S.
      • Mao J.-H.
      • et al.
      An interferon signature identified by RNA-sequencing of mammary tissues varies across the estrous cycle and is predictive of metastasis-free survival.
      Advancing age and whether a mouse has experienced pregnancy can modify relative epithelial cell content.
      • Raafat A.
      • Strizzi L.
      • Lashin K.
      • et al.
      Effects of age and parity on mammary gland lesions and progenitor cells in the FVB/N-RC mice.
      Mammary gland–targeted hormonal interventions, such as overexpression of estrogen receptor α or aromatase in mammary epithelial cells, induce pathologic changes, including increased lobular density, hyperplastic alveolar nodules, and cancer.
      • Dabydeen S.A.
      • Kang K.
      • Díaz-Cruz E.S.
      • Alamri A.
      • Axelrod M.L.
      • Bouker K.B.
      • Al-Kharboosh R.
      • Clarke R.
      • Hennighausen L.
      • Furth P.A.
      Comparison of tamoxifen and letrozole response in mammary preneoplasia of ER and aromatase overexpressing mice defines an immune-associated gene signature linked to tamoxifen resistance.
      ,
      • Díaz-Cruz E.S.
      • Sugimoto Y.
      • Gallicano G.I.
      • Brueggemeier R.W.
      • Furth P.A.
      Comparison of increased aromatase versus ERα in the generation of mammary hyperplasia and cancer.
      ,
      • Frech M.S.
      • Halama E.D.
      • Tilli M.T.
      • et al.
      Deregulated estrogen receptor alpha expression in mammary epithelial cells of transgenic mice results in the development of ductal carcinoma in situ.
      This study evaluates a solution that fills the need for mouse mammary gland density assessment across age, reproductive senescence, and genetically or chemically induced pathophysiology. It can allow direct comparison of dense lobular morphologies to morphologies with clear branching patterns. Although sophisticated methods exist for the examination of density data, the process used in this study provides a streamlined and efficient method to replace laborious and subjective visual density scoring while still relying on standard imaging equipment and techniques. Screening of mouse mammary glands for density can be used to identify mice at higher risk for cancer, just as mammography is used for women. Having a process to accomplish this for mice can aid an investigator in identifying which glands in an experiment should be targeted for focused in-depth study of cancer pathophysiology. Because density screening can facilitate identification of glands with nonpalpable cancer, the process enables study of smaller cancers at earlier stages in their generation.
      Adequate denoising is required for image analysis.
      • Fan L.
      • Zhang F.
      • Fan H.
      • Zhang C.
      Brief review of image denoising techniques.
      In the case of normal mammary gland whole mounts, higher-density branching and lobular growth patterns lie within a low-density fat pad. Higher-density pathologic findings include hyperplastic alveolar nodules, cancers, and a fibrotic and/or proliferative fat pad.
      • Frech M.S.
      • Torre K.M.
      • Robinson G.W.
      • Furth P.A.
      Loss of cyclin D1 in concert with deregulated estrogen receptor alpha expression induces DNA damage response activation and interrupts mammary gland morphogenesis.
      Gaussian denoising is one of the most commonly used methods for image denoising and results in the removal of typical speckle noises frequently present in images.
      • Kumar N.
      • Nachamai M.
      Noise removal and filtering techniques used in medical images.
      Ridge operators are used in image analysis to detect ridge-like structures, such as ducts.
      • Meijering E.
      • Jacob M.
      • Sarria J-Cf
      • Steiner P.
      • Hirling H.
      • Unser M.
      Design and validation of a tool for neurite tracing and analysis in fluorescence microscopy images.
      This study evaluates a bifurcated approach using two different image processing approaches, gaussian denoising and application of a ridge operator, for the analysis of mouse mammary gland density using digital images of mammary gland whole mounts.

      Materials and Methods

      Mouse Models

      The animal research protocol was approved by the Georgetown University Institutional Animal Care and Use Committee. All regulations concerning the use of animals in research were adhered to. To ensure that a range of mammary gland whole mount morphologies were available for testing of the digital imaging protocol, whole mount images from C57Bl/6 genetically engineered mice over a broad age and weight range under different transgene induction conditions were used. Two different genetically engineered mouse models were used, one with conditional overexpression of Esr1
      • Frech M.S.
      • Halama E.D.
      • Tilli M.T.
      • et al.
      Deregulated estrogen receptor alpha expression in mammary epithelial cells of transgenic mice results in the development of ductal carcinoma in situ.
      and one with conditional overexpression of CYP19A1.
      • Díaz-Cruz E.S.
      • Sugimoto Y.
      • Gallicano G.I.
      • Brueggemeier R.W.
      • Furth P.A.
      Comparison of increased aromatase versus ERα in the generation of mammary hyperplasia and cancer.
      All whole mounts generated from mice that did not reach target age end points of a parent study investigating the impact of age on mammary gland disease and antihormonal response were included in this image analysis study. Only whole mounts that were improperly stained or incorrectly dissected with skin or muscle tissue inclusion or other artifact in the imaging field were excluded from analysis. The numbers of whole mounts from each genetically engineered mouse model and transgene induction condition included were no transgene induction (Esr1 n = 18, CYP19A1 n = 8), transgene induction at 12 months of age (Esr1 n = 43, CYP19A1 n = 59), and transgene induction at 18 months of age (Esr1 n = 2, CYP19A1 n = 4). Reasons for necropsy of the mice included planned end points for this density study at 4, 5, and 27 months of age (n = 20) and necropsies for humane reasons per approved animal protocol between 9 and 29 months of age (n = 114). The total number of whole mounts in the study (n = 134), the means ± SE mouse age (19 ± 0.6 months), and the means ± SEM mouse weight (29 ± 0.8 kg) were determined. Sample sizes were determined statistically before experimentation for analysis of correlations with visual density ranking and specific mammary gland morphologies. Samples were assigned unique numerical identifiers and processed in order of availability for the study, independent of genotype and cohort identity, and data were acquired blindly and recorded by the unique numerical identifier.

      Mammary Gland Histologic Evaluation

      Carmine alum–stained inguinal mammary gland whole mounts were prepared according to standard procedures.
      • Tolg C.
      • Cowman M.
      • Turley E.A.
      Mouse mammary gland whole mount preparation and analysis.
      Whole mounts were visualized using a Nikon Eclipse E800M microscope (Nikon Corporation, Melville, NY) at 0.5×, 1×, and 4× to evaluate secondary versus tertiary branching and presence of lobular growth and hyperplastic alveolar nodules as previously described.
      • Alothman S.J.
      • Wang W.
      • Goerlitz D.S.
      • Islam M.
      • Zhong X.
      • Kishore A.
      • Azhar R.I.
      • Kallakury B.V.
      • Furth P.A.
      Responsiveness of Brca1 and Trp53 deficiency–induced mammary preneoplasia to selective estrogen modulators versus an aromatase inhibitor in Mus musculus.
      ,
      • Nakles R.E.
      • Kallakury B.V.S.
      • Furth P.A.
      The PPARγ agonist efatutazone increases the spectrum of well-differentiated mammary cancer subtypes initiated by loss of full-length BRCA1 in association with TP53 haploinsufficiency.
      Two to four blinded observers evaluated each parameter independently (B.L.R., V.M., W.W., and/or P.A.F.). When observers reported different findings, the majority finding was adopted for the final reading. When there was a tie in readings, a secondary reading was performed for final assessment (P.A.F.) Mammary gland histology was read blindly on hematoxylin and eosin–stained sections of formalin-fixed, paraffin-embedded tissue with scoring for the presence or absence of invasive cancer and/or a proliferative or fibrotic stroma using a Nikon Eclipse E800M microscope at ×40 (P.A.F.). Evaluation confirmed that the population of whole mounts included in the image analysis study showed a range of normal and pathological findings: lobular and/or lobuloalveolar (lobular) growth (n = 75), no lobular growth (n = 59), secondary branching (n = 55), secondary and tertiary (tertiary) branching (n = 79), hyperplastic alveolar nodule (HAN) present (n = 46), no HANs (n = 88), mammary cancer present (n = 34), no mammary cancer (n = 105), normal fat pad (n = 106), and fibrotic or proliferative fat pad (n = 28).

      Whole Mount Digital Image Preparation

      High-resolution red, green, blue (RGB) TIF images (14.1 MB, 2560 × 1920 pixels) were obtained using a Nikon Eclipse E800M microscope equipped with a DMX1200 camera with the NIS Elements imaging software version 4.30 (Nikon Corporation, Tokyo, Japan) using a 0.5× objective. Inguinal mammary gland whole mount structures distal to the lymph node were uniformly imaged. When this area could not be imaged within one field, two or more overlapping images were taken digitally and combined into a single image using Adobe Photoshop 2021 version 22.3.0 (Adobe Inc., San Jose, CA). Digital imaging requirements included a uniformly well-focused image with an absence of muscle or skin fascia in the central region of the whole mount, staining, or mounting artifacts, all of which compromised the density measurement. Once an RGB digital image was obtained, images were manually cropped such that the top and bottom borders of the image aligned with the innermost border of the fat pad generating representative rectangular images that lacked any surrounding white space.

      Program

      A Python program was written to use scikit-image version 0.16.2 to grayscale RGB digital images of mouse mammary gland whole mounts.
      • van der Walt S.
      • Schönberger J.L.
      • Nunez-Iglesias J.
      • Boulogne F.
      • Warner J.D.
      • Yager N.
      • Gouillart E.
      • Yu T.
      The scikit-image Contributors
      scikit-image: image processing in Python.
      After conversion to gray scale, the program automatically cropped 30% off the distal side of the image and 10% off the proximal side of the image to remove regions frequently confounded by artifacts, specifically skeletal muscle on the distal end and dark staining around the lymph node, to center measurement of relative mammary gland density on the midgland distal to the lymph node. Initial mean pixel intensity was measured with NumPy version 1.19.2.
      • Harris C.R.
      • Millman K.J.
      • van der Walt S.J.
      • et al.
      Array programming with NumPy.
      Intensity values range from 0 (all black image) to 255 (all white image). Images with mean pixel intensity <89.99 were sorted to a gaussian denoising method,
      • van der Walt S.
      • Schönberger J.L.
      • Nunez-Iglesias J.
      • Boulogne F.
      • Warner J.D.
      • Yager N.
      • Gouillart E.
      • Yu T.
      The scikit-image Contributors
      scikit-image: image processing in Python.
      and images with mean pixel intensity >89.99 were sorted to a ridge operator.
      • Meijering E.
      • Jacob M.
      • Sarria J-Cf
      • Steiner P.
      • Hirling H.
      • Unser M.
      Design and validation of a tool for neurite tracing and analysis in fluorescence microscopy images.
      Gaussian denoising parameters were optimized using three images visually identified as low, medium, and high density. These images were processed through the denoising portion of the assay. Each parameter of the bilateral denoising, color, and spatial closeness has a programmable SD value in which larger values result in pixels with greater differences in size or color value being averaged. Beginning with the color parameter, programmable as a float value of 0.0 to 1.0 (default is none), images were tested at a range of values, whereas the spatial parameter was left by default at 15, and visually analyzed for the optimal parameter for removing background features, such as staining of the fat pad. The same analysis was conducted for the spatial parameter, programmable as positive values (default was 1), with the visually optimized color parameter. The optimal pairing was identified as a color parameter of 0.1 and a spatial parameter of 25. A method for neurite tracing was applied and found to accurately identify branching structures in the whole mounts.
      • Meijering E.
      • Jacob M.
      • Sarria J-Cf
      • Steiner P.
      • Hirling H.
      • Unser M.
      Design and validation of a tool for neurite tracing and analysis in fluorescence microscopy images.
      This process was visually optimized by applying morphology closing and removing small object functions from the scikit image to remove nonductal structures from the background of the whole mount before application of the ridge operator.
      • van der Walt S.
      • Schönberger J.L.
      • Nunez-Iglesias J.
      • Boulogne F.
      • Warner J.D.
      • Yager N.
      • Gouillart E.
      • Yu T.
      The scikit-image Contributors
      scikit-image: image processing in Python.
      After denoising, the program calculated arithmetic mean intensity using NumPy version 1.19.2
      • Harris C.R.
      • Millman K.J.
      • van der Walt S.J.
      • et al.
      Array programming with NumPy.
      as a proxy for gland density in which epithelium and dense stroma appeared darker and free space and fat pad appeared lighter. A quality control check was performed by manual visual inspection of images to confirm that inadvertent background had not been included or ductal structures excluded. Images that were determined to have been improperly sorted into the gaussian denoising pathway but that had clear linear ductal structures were reentered into the program with the denoising pathway disabled, ensuring processing via ridge operator. With use of the initial mean pixel density sort set at 89.99, all the images optimally processed by gaussian denoising were correctly sorted, but approximately 20% of the images with initial mean pixel intensity <89.99 required a ridge operator for optimal denoising. The 89.99 threshold was chosen after visual optimization and mean pixel intensity measurement. A subset of images was subjected to both processing methods, and their preprocessing gray scale mean pixel intensities were recorded. The output images were all checked for quality control to verify which path gave the most accurate representation of the ductal structure. All images that were correctly processed by the ridge filter had a mean pixel intensity of >90.0 (minimum, 90.3), and most images correctly processed by gaussian denoising had initial mean pixel intensity values <90.0. Accordingly, 89.99 was chosen so that corrective sorting only needed to occur in one direction, moving images from the gaussian denoising pathway to the ridge filter pathway and not vice versa to minimize operator effort and bias. Manual cropping of red, green, blue images of carmine alum–stained mammary gland whole mount images was followed by automated cropping within the program to generate rectangular images of each whole mount that excluded lymph node and irregular edge structures. The program then converted the image to grayscale and performed an initial measurement of mean pixel intensity. When mean pixel intensity was ≤89.99, images were automatically sorted to gaussian denoising for isolation of ductal structure from background, and a final mean pixel intensity (density score) was determined. This process was followed by a manual quality control step. If an image was determined to have a predominantly linear rather than spherical structure format at this step, it was manually re-sorted to application of a ridge operator for image processing. When mean pixel intensity was >89.99, images were automatically sorted to ridge operator image processing, followed by a manual quality control step. The program includes generation of histograms for the processed images, created using Matplotlib version 3.4.1.
      • Hunter J.D.
      Matplotlib: a 2D graphics environment.
      In addition to the arithmetic mean of pixel intensities, the program also reports median, using Numpy version 1.19.2, and geometric mean, using SciPy version 1.8.1, for processed images.
      • Harris C.R.
      • Millman K.J.
      • van der Walt S.J.
      • et al.
      Array programming with NumPy.
      ,
      • Virtanen P.
      • Gommers R.
      • Oliphant T.E.
      • et al.
      SciPy 1.0 Contributors
      SciPy 1.0: fundamental algorithms for scientific computing in Python.
      For data presented here, arithmetic mean was used after testing whether statistically significant associations were altered by selecting one or another of these output measures. No difference in the comparisons presented here was found dependent on output score type, that is, arithmetic mean, median mean, or geometric mean. However, this functionality was included in the program because it is possible that future studies would preferentially need one or the other of these output measures. The program code is presented in Supplemental Code S1.

      Comparison of Relative Mammary Gland Density Measurement with Manual Visual Relative Mammary Gland Density Assessment

      To determine whether the relative mammary gland density readings determined through the two different pathways correlated with visual density rankings, RGB whole mount digital images from mice with transgene induction at 12 months of age were selected from each genotype (Esr1 n = 29; CYP19A1 n = 33) for manual visual ranking by two blinded independent observers (B.L.R. and P.A.F.). Each observer ranked the images from least to most dense (1 to 29 for Esr1 and 1 to 33 for CYP19A1). A summative visual ranking score was obtained by combining the ranking number of each independent investigator.

      Comparison of Relative Mammary Gland Density Measurement with Histologic Findings

      To determine how relative mammary gland density readings correlated with histologic findings, relative density scores were compared between whole mounts with secondary versus tertiary branching and presence versus absence of lobular growth, HANs, invasive cancer, and abnormal fibrotic or proliferative stroma.

      Statistical Analysis

      Linear regression analyses were used to compare relative density measurements with manual visual density assessment. U-tests were used to compare relative density measurements in whole mounts with secondary versus secondary and tertiary (tertiary) branching and the absence versus presence of lobular growth, fibrotic or proliferative fat pad, mammary cancer and hyperplastic alveolar nodules. P < 0.05 was considered statistically significant (GraphPad Prism version 9.3.1, GraphPad Software Inc., San Diego, CA).

      Results

      Sorting of Mammary Gland Whole Mount Images into Gaussian Denoising or Ridge Operator Processing Pathways

      The bifurcated program sorted whole mount morphologies (n = 134) into two different image processing approaches using an initial mean pixel intensity measurement following manual and automated image cropping and grayscale conversion (Figure 1A). Application of an initial 89.99-pixel threshold enabled automated sorting into the two different pathways. All images optimally processed using gaussian denoising were sorted correctly using this threshold. Approximately 20% of images with an initial mean pixel intensity <89.99 were assessed at the quality control step to need ridge operator processing and required manual sorting to the ridge operator sequence. The manual visual quality control step assessed whether nonductal structure background was included in the processed image. When the predominant ductal structures included more spherical lobular growth, these images could be processed through gaussian denoising without background contamination (Figure 1B and Supplemental Figure S1A). However, if significant background was present in a whole mount image that exhibited a more linear growth pattern, initial mean pixel intensity could decrease to <89.99 with inclusion of this background in the relative density score if processed by gaussian denoising (Figure 1C). In contrast, the same image processed using a ridge operator would not include this artifactual background (Figure 1C). Most linear-structured glands without background contamination yielded initial mean pixel intensity scores >89.99 and were correctly sorted into ridge operator processing using the automated decision step (Figure 1D). Of the 134 images that were assessed using the bifurcated program, 57 were optimally processed using gaussian denoising and 77 were optimally processed using a ridge operator. Arithmetic mean density scores were derived from the relative pixel count (Supplemental Figure S1B). Histograms show that the final appropriately processed images produce narrow pixel distributions (Figure 1, B–D). Images that exhibited initial mean pixel intensity measurements ≤89.99 but because of image characteristics were more optimally processed through a ridge filter demonstrated a broader pixel distribution with gaussian denoising (Figure 1C).
      Figure thumbnail gr1
      Figure 1Density measurement assessment using gaussian denoising or a ridge operator for image processing. A: Flowchart of digital image processing and density measurement program. Directional arrows indicate passage sequence through the different steps of the digital image processing and density measurement program. Red slanted arrows indicate automated decision points, and the orange slanted arrow indicates a manual decision point. B: Representative image processing of a carmine alum–stained mammary gland whole mount image through gaussian denoising. Histogram of image processed by gaussian denoising shown. C: Representative image processing of a carmine alum–stained mammary gland whole mount image initially processed through gaussian denoising with re-sorting to a ridge operator for optimal density measurement after the quality control step. White arrow indicates darker artifactual background strip that was included in pixel count when gaussian denoising is applied. Histograms of images processed by gaussian denoising or by ridge operator shown. D: Representative image processing of a carmine alum–stained mammary gland whole mount image through a ridge operator. Images of carmine alum–stained mammary gland whole mounts were taken with a Nikon Eclipse microscope. Histogram of image processed by ridge operator shown. All whole mount images are to the same scale. Scale bar = 1000 μm (B). Original magnification, ×0.5.

      Relative Mean Density Scores Rendered by the Bifurcated Program across a Range of Mouse Ages and Physiologic States

      The program was used to analyze whole mount images from mice 4 to 29 months of age and assigned arithmetic mean density scores ranging from 14 (highest density) to 220 (lowest density). Simple secondary branching characteristic of estrus had relatively lower density (score 165) (Figure 2A), whereas lobuloalveolar growth found during diestrus had high density (score 36) (Figure 2B). Mammary gland ductal changes are a continuum through the estrus cycle. Thinner secondary branching structures rendered a lower density score (score 220) (Figure 2C), whereas the presence of both secondary and tertiary branching increased the relative density (score 166) (Figure 2D). Density measurements enabled assessment of the range of density found in aging mice with disease-inducing transgene expression from a relative mean density score of 35 when dense lobuloalveolar growth was present (Figure 2E) to a relative mean density score of 156 in a gland demonstrating only secondary and tertiary branching (Figure 2F). Histograms demonstrate differences in distribution of pixel intensities among the different processed images.
      Figure thumbnail gr2
      Figure 2Comparison of mammary gland structures present in images processed by gaussian denoising or a ridge operator. A: Ridge operator processing of a carmine alum–stained mammary gland whole mount image with predominantly secondary branching from a 4-month–old mouse. Relative density score is 165. B: Gaussian denoising processing of a carmine alum–stained mammary gland whole mount image with predominantly lobuloalveolar growth from a 9-month–old mouse. Relative density score is 36. C: Ridge operator processing of a carmine alum–stained mammary gland whole mount image with predominantly secondary branching from a 14-month–old mouse. Relative density score is 220. D: Ridge operator processing of a carmine alum–stained mammary gland whole mount image with both secondary and tertiary branching from a 20-month–old mouse. Relative density score is 166. E: Gaussian denoising processing of a carmine alum–stained mammary gland whole mount image with predominantly lobuloalveolar growth from a 22-month–old mouse. Relative density score is 35. F: Ridge operator processing of a carmine alum–stained mammary gland whole mount image with both secondary and tertiary branching from a 25-month–old mouse. Relative density score is 156. Corresponding regions of original red, green, blue (RGB) image and ridge operator output of grayscale image are shown in panels A–F. Histograms of images processed by ridge operator (A, C, D, and F) or gaussian denoising (B and E). Images of carmine alum–stained mammary gland whole mounts taken with a Nikon Eclipse microscope. All whole mount images are to the same scale. Scale bar = 1000 μm (A). Original magnification, ×0.5.

      Correlation of the Relative Mean Density Scores Rendered by the Bifurcated Program with Manual Visual Density Assessment

      The two different genetically engineered mouse models (Esr1 and CYP19A1) used for development of the density assessment program were analyzed separately for comparison of digital and visual densities. For the assessment presented here, only Esr1 (n = 29) and CYP19A1 (n = 33) mice with transgene induction at 12 months of age were used. Blinded visual rankings were compared with the digital density results and analyzed by simple linear regression (GraphPad Prism version 9.3.1). Although the calculated slope of Esr1 (Figure 3A) differed slightly from that of CYP19A1 (Figure 3B) mice, significant correlations between digital and visual density rankings were obtained for both models (P < 0.0001). In both models, lower density was associated with predominantly secondary branching (Figure 3, C–F), appearance of tertiary and secondary branching with slightly higher density (Figure 3, G–H), and lobuloalveolar growth with highest density (Figure 3, I–L). Differences in distribution of pixel intensities between processed images are shown in histograms for each image.
      Figure thumbnail gr3
      Figure 3Correlation of digital density measurements with manual visual density rankings. Linear regression plot comparing relative density scores with manual visual rankings of carmine alum–stained mammary gland whole mount images for Esr1 mice with transgene induction at 12 months of age (A) and CYP19A1 mice with transgene induction at 12 months of age (B). Corresponding regions of original red, green, blue (RGB) image and ridge operator output from a carmine alum–stained mammary gland whole mount image from a 28-month–old Esr1 mouse with a density score of 212 (C); a 17-month–old CYP19A1 mouse with a density score of 192 (D); a 28-month–old Esr1 mouse with a relative density score of 193 (E); a 26-month–old CYP19A1 mouse with a density score of 182 (F); a 13-month–old Esr1 mouse with a density score of 147 (G); a 21-month–old CYP19A1 mouse with a density score of 151 (H); a 27-month–old Esr1 mouse with a density score of 59 (I); a 27-month–old CYP19A1 mouse with a density score of 27 (J); a 17-month–old Esr1 mouse with a density score of 47 (K); and a 25-month–old CYP19A1 mouse with a density score of 21 (L). Histograms of images processed by ridge operator (CH, lower panels) or gaussian denoising (IL, lower panels) shown below processed image. Images of carmine alum–stained mammary gland whole mounts taken with a a Nikon Eclipse microscope. Linear regression analysis Esr1: goodness of fit: R2 = 0.8864, Sy.x = 5.697. F = 210.6 DFn, DFd = 1.27, P < 0.0001. Equation: Y = −0.2616∗X + 60.42. n = 29. Linear regression analysis CYP19A1: goodness of fit: R2 = 0.8314, Sy.x = 7.822. F = 152.8 DFn, DFd = 1.31, P < 0.0001. Equation: Y = −0.3131 ∗ X + 73.91. n = 33. All whole mount images are to the same scale. Scale bar = 1000 μm (C). Original magnification, ×0.5.

      Significant Associations with Tertiary Branching, Lobular Growth, Fibrotic Stroma, and Presence of Mammary Cancer

      To further evaluate correlations between density measurements and specific features of the individual whole mounts, digital density measurements of all 134 whole mounts in the study were analyzed for correlation with branching patterns and the presence or absence of lobular growth, a fibrotic fat pad, mammary cancer, or HANs (U-test, two-tailed, GraphPad Prism version 9.3.1). Presence of only secondary branching was significantly associated with lower density than the combination of secondary and tertiary branching as indicated by the higher mean density scores (P < 0.0001) (Figure 4A). Similarly, lack of lobular growth was significantly associated with lower and presence of lobular growth with higher density (P < 0.0001) (Figure 4B). The presence of a fibrotic fat pad with proliferative stroma was significantly associated with lower density (P < 0.001) (Figure 4C). When mammary cancer was present, it was significantly associated with higher density (P < 0.05) (Figure 4D). Of interest, detection of HANs was significantly correlated with lower density, perhaps because they could be more readily appreciated in a lower density gland (P < 0.01) (Figure 4E).
      Figure thumbnail gr4
      Figure 4Higher density correlates with the presence of lobular growth, tertiary branching, fibrotic mammary fat pad, and mammary cancer. A: Scatterplot illustrating relative density scores in whole mounts with and without lobular growth. Means ± SEM: of 150 ± 6 indicate no lobular growth. Means ± SEM of 86 ± 6 indicate lobular growth. Carmine alum–stained mammary gland whole mount images of secondary and tertiary with secondary ductal branching are shown. B: Scatterplot illustrating relative density scores in whole mounts with secondary compared with tertiary branching. Means ± SEM of 145 ± 7 indicate secondary branching. Means ± SEM of 93 ± 6 indicate tertiary branching. Carmine alum–stained mammary gland whole mount images of ductal structures without and with lobular growth are shown. C: Scatterplot illustrating relative density scores in whole mounts with normal mammary fat pads compared with subnormal fibrotic or proliferative fat pads. Means ± SEM of 124 ± 5 indicate mammary fat pad. Means ± SEM of 79 ± 11 indicate fibrotic proliferative fat pad. Carmine alum–stained mammary gland whole mount images of normal and abnormal fibrotic or proliferative fat pad are shown. D: Scatterplot illustrating relative density scores in whole mounts without and with mammary cancer. Means ± SEM of 121 ± 7 indicate no mammary cancer. Means ± SEM of 87 ± 92 indicate mammary cancer. Carmine alum–stained mammary gland whole mount images of mammary glands without and with mammary cancer are shown. E: Scatterplot illustrating relative density scores in whole mounts with and without hyperplastic alveolar nodule (HANs). Means ± SEM of 104 ± 6 indicates no HANs. Means ± SEM of 133 ± 8 indicates HANs. Carmine alum–stained mammary gland whole mount images of mammary glands with and without a HAN are shown. Images of carmine alum–stained mammary gland whole mounts taken with a Nikon Eclipse microscope. n = 134 all comparisons. ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001, and ∗∗∗∗P < 0.0001. All whole mount images are to the same scale. Scale bar = 1000 μm (A). Original magnification, ×0.5.

      Discussion

      This study demonstrated the utility of a bifurcated approach for the assessment of mammary gland density using images of mammary gland whole mounts. A single image processing approach was not able to execute optimal denoising across different mammary gland morphologies. For example, although a ridge operator adequately handled whole mount images with simpler linear branching patterns, it failed to isolate epithelial growth away from background when dense lobular growth and extensive tertiary branching were present. For these glands, application of gaussian denoising was the optimal image processing approach for epithelial density measurement. To optimize a digital flow of image analysis, it was determined that an initial read of mean pixel intensity could be used for initial sorting into the two image processing pathways. Execution of this approach enabled a quantitative measurement of mammary gland whole mount density over a wide range of normal and abnormal findings. Higher densities were associated with lobular growth, tertiary branching, fibrotic and proliferative stroma, and cancers. Inclusion of code that enabled histogram generation for each processed image enabled comparisons of pixel intensity distributions across different images and the two different processing approaches.
      Similar to measurements of mammary gland density in women, positive associations were found between higher density and the presence of mammary neoplasia or fibrosis.
      • Sherratt M.J.
      • McConnell J.C.
      • Streuli C.H.
      Raised mammographic density: causative mechanisms and biological consequences.
      HAN detection was correlated with relatively lower mammary gland density measurements. This observation is reminiscent of the problem in assessing the presence of cancer in women with higher mammary density, as observed in younger women when masking of tumors in the setting of high mammary gland density occurs.
      • Nazari S.S.
      • Mukherjee P.
      An overview of mammographic density and its association with breast cancer.
      ,
      • Mainprize J.G.
      • Alonzo-Proulx O.
      • Alshafeiy T.I.
      • Patrie J.T.
      • Harvey J.A.
      • Yaffe M.J.
      Prediction of cancer masking in screening mammography using density and textural features.
      HANs are conventionally defined as single discrete lesions distinct from normal ductal and lobular structures. Some of the glands imaged in this study showed high density that correlated with nearly uniform lobular growth or contained cancers that masked identification of discrete HANs.
      Mammary density is an important risk factor for human cancer development, yet it is not routinely included as a metric in rodent mammary gland whole mount analysis because of challenges in quantification. The method analyzed here can contribute to increasing the availability of mammary gland density data in future studies, complementing related but different metrics to measure branching density
      • Fernandez-Gonzalez R.
      • Barcellos-Hoff M.H.
      • Ortiz-de-Solórzano C.
      Quantitative image analysis in mammary gland biology.
      • Stanko J.P.
      • Easterling M.R.
      • Fenton S.E.
      Application of Sholl analysis to quantify changes in growth and development in rat mammary gland whole mounts.
      • McGinley J.N.
      • Thompson H.J.
      Quantitative assessment of mammary gland density in rodents using digital image analysis.
      and in so doing promoting better understanding of the role of mammary gland density in cancer. This method has the added benefit of being applicable across a wide range of ages, morphologies, and disease, enabling its use even where clear branching patterns are not present. Previously developed imaging protocols have centered on evaluation of younger mice. The data presented here document development of an imaging protocol that can be successfully applied to aging mice up to 29 months of age. The program includes modifiable parameters, such as the sorting threshold and denoising strength, promoting broad applicability for routine integration into mouse model studies wherever researchers have access to standard two-dimensional whole mount images. Future efforts to automate the preprocessing steps used here would further streamline the process.
      In this study, increased mammary density was correlated with the presence of mammary cancer. In future projects, one use of density measurement would be to incorporate it into the experimental flow of a study on mammary cancer generation. It could enable sorting of samples into higher and lower probabilities of cancer, which could be beneficial from a practical cost and effort perspective. For example, instead of sectioning all glands for cancer screening, density screening could be used to limit that step only to the higher-risk samples. Similarly, it could assist in select samples for primary cell culture. There are means to perform initial processing of tissue for histology or cell culture that could then be held in abeyance for final processing until after mammary density could be evaluated. This approach would allow studies to proceed more quickly and cost-effectively. Obtaining relative density scores enabled statistical comparisons with different normal and abnormal growth parameters. Because scores were obtained from a digital program, intraobserver and interobserver errors were eliminated, which is also one of the intended benefits of imaging analysis software for human mammogram assessment.
      • Kerlikowske K.
      • Scott C.G.
      • Mahmoudzadeh A.P.
      • et al.
      Automated and clinical breast imaging reporting and data system density measures predict risk for screen-detected and interval cancers: a case-control study.
      Results from this study showed how automated assessment of breast density complements visual whole mount assessment, just as advocated for human mammography.
      The pairwise comparisons between groups revealed significant associations between density and specific physiologic and pathophysiologic features of the mammary gland whole mounts. They also illustrated the bimodal distribution of the data, with a split in density scores occurring along the bifurcated method. If desired, application of a log transformation of the data could be used for normalization of data derived from the density assessment program. Here, even with a bimodal distribution, there was high correlation between relative density measurements and visual assessment, and statistically significant differences between pairwise comparisons were not altered by data normalization. Because the data distribution was nonnormal, nonparametric tests were used for analysis of differences between groups.
      To implement the bifurcated processing strategy, a differentiating point had to be decided. The process described in this study sought to maximize automation by using a numerical decision point rather than a visual one, which was accomplished with the 89.99 cutoff value for mean pixel intensity of the prefiltered gray images. The number 89.9 was chosen as the threshold for this decision point with the knowledge that it would sort some images into gaussian denoising that would be optimally handled by a ridge operator but would correctly sort all images optimally handled by gaussian denoising into the correct pathway. The application of a numerical cutoff rather than visually sorting was the most efficient approach because there was no invariable correlation between a specific structural element and the optimal imaging processing pathway.
      The program described here used arithmetic mean for the cutoff but could be easily modified to use median or geometric mean. A different approach for sorting would be to discriminate the use of histograms of the processed images. As shown here, when an image with a score <89.99 yields a false density score because of incorrect background inclusion, the histogram shows a broader pixel intensity distribution compared with when it is correctly processed using ridge filtering. The program was designed so that individual operators could optimize the program to their own data set. For example, use of a different set of whole mount images with a different distribution of biological characteristics might require a different cutoff than the 89.99 setpoint used here. Use of median or geometric mean as the sorting measure would require recalibration of the numerical cutoff for optimal program functioning.
      The code used for the approach was developed from publicly available sources using previously validated methods for image processing. The flow used here can be implemented without the need for specialized equipment or additional software purchase and is designed as an efficient replacement for visual density scoring. Development of this bifurcated flow process enabled density measurements across a wide range of whole mount findings, both normal and abnormal.
      Limitations of the program include the requirement for meticulous whole mount preparation to avoid artifacts that could not be corrected with these image processing techniques and operator involvement to perform manual quality check after machine sorting by initial mean pixel intensity to enable optimal image processing for each image.

      Acknowledgment

      We thank Zachary Susswein for his valuable help and advice on the Python code.

      Author Contributions

      B.L.R. conceptualized the study; B.L.R. and P.A.F. designed the study; B.L.R., B.P.R., V.M., W.W., and P.A.F. acquired the data; B.L.R. and P.A.F. analyzed and interpreted the data and prepared the initial manuscript; B.L.R., B.P.R., V.M., W.W., and P.A.F. reviewed the manuscript before submission and acceptance of authorship; P.A.F. secured funding. P.A.F. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

      Supplemental Data

      Figure thumbnail figs1
      Supplemental Figure S1A: Illustration of a mammary gland whole mount image exhibiting spherical lobular growth that would be optimally processed using gaussian denoising demonstrating how processing through a ridge operator would yield artifactual density scoring. Original gray, cropped, gaussian denoising and ridge operator processed images are shown. When the ridge filter was applied, it failed to find the branching type structures it was designed to isolate. This failure produced an erroneous output image with large white spaces that would result in an inaccurately high (less dense) density score when the score more accurately would be lower (more dense). B: Scatterplot illustrating gaussian denoising arithmetic mean density scores compared with ridge operator arithmetic mean density scores for mammary gland whole mounts analyzed in the study. Blue dots are the gray cropped images correctly sorted to gaussian denoising pathway using the 89.99-pixel threshold. Pink dots are the gray cropped images correctly sorted to the ridge operator pathway using the 89.99-pixel threshold. Orange dots are the gray cropped images incorrectly sorted to gaussian denoising pathway using the 89.99-pixel threshold that were manually corrected to the ridge operator pathway for optimal processing.

      References

        • Fuller M.S.
        • Lee C.I.
        • Elmore J.G.
        Breast cancer screening: an evidence-based update.
        Med Clin North Am. 2015; 99: 451-468
        • Burton A.
        • Maskarinec G.
        • Perez-Gomez B.
        • et al.
        Mammographic density and ageing: a collaborative pooled analysis of cross-sectional data from 22 countries worldwide.
        PLoS Med. 2017; 14: e1002335
        • Carney P.A.
        • Miglioretti D.L.
        • Yankaskas B.C.
        • et al.
        Individual and combined effects of age, breast density, and hormone replacement therapy use on the accuracy of screening mammography.
        Ann Intern Med. 2003; 138: 168-175
        • Mendes J.
        • Matela N.
        Breast cancer risk assessment: a review on mammography-based approaches.
        J Imaging. 2021; 7: 98
        • Boyd N.F.
        • Martin L.J.
        • Yaffe M.J.
        • Minkin S.
        Mammographic density and breast cancer risk: current understanding and future prospects.
        Breast Cancer Res. 2011; 13: 223
        • Jakes R.
        • Duffy S.
        • Ng F.
        • Gao F.
        • Ng E.
        Mammographic parenchymal patterns and risk of breast cancer at and after a prevalence screen in Singaporean women.
        Int J Epidemiol. 2000; 29: 11-19
        • Huo Z.
        • Giger M.L.
        • Wolverton D.E.
        • Zhong W.
        • Cumming S.
        • Olopade O.I.
        Computerized analysis of mammographic parenchymal patterns for breast cancer risk assessment: feature selection.
        Med Phys. 2000; 27: 4-12
        • Dabydeen S.A.
        • Kang K.
        • Díaz-Cruz E.S.
        • Alamri A.
        • Axelrod M.L.
        • Bouker K.B.
        • Al-Kharboosh R.
        • Clarke R.
        • Hennighausen L.
        • Furth P.A.
        Comparison of tamoxifen and letrozole response in mammary preneoplasia of ER and aromatase overexpressing mice defines an immune-associated gene signature linked to tamoxifen resistance.
        Carcinogenesis. 2015; 36: 122-132
        • Alothman S.J.
        • Wang W.
        • Goerlitz D.S.
        • Islam M.
        • Zhong X.
        • Kishore A.
        • Azhar R.I.
        • Kallakury B.V.
        • Furth P.A.
        Responsiveness of Brca1 and Trp53 deficiency–induced mammary preneoplasia to selective estrogen modulators versus an aromatase inhibitor in Mus musculus.
        Cancer Prev Res. 2017; 10: 244-254
        • Nakles R.E.
        • Kallakury B.V.S.
        • Furth P.A.
        The PPARγ agonist efatutazone increases the spectrum of well-differentiated mammary cancer subtypes initiated by loss of full-length BRCA1 in association with TP53 haploinsufficiency.
        Am J Pathol. 2013; 182: 1976-1985
        • Fernandez-Gonzalez R.
        • Barcellos-Hoff M.H.
        • Ortiz-de-Solórzano C.
        Quantitative image analysis in mammary gland biology.
        J Mammary Gland Biol Neoplasia. 2004; 9: 343-359
        • Stanko J.P.
        • Easterling M.R.
        • Fenton S.E.
        Application of Sholl analysis to quantify changes in growth and development in rat mammary gland whole mounts.
        Reprod Toxicol. 2015; 54: 129-135
        • McGinley J.N.
        • Thompson H.J.
        Quantitative assessment of mammary gland density in rodents using digital image analysis.
        Biol Proced Online. 2011; 13: 4
        • Snijders A.M.
        • Langley S.
        • Mao J.-H.
        • et al.
        An interferon signature identified by RNA-sequencing of mammary tissues varies across the estrous cycle and is predictive of metastasis-free survival.
        Oncotarget. 2014; 5: 4011-4025
        • Raafat A.
        • Strizzi L.
        • Lashin K.
        • et al.
        Effects of age and parity on mammary gland lesions and progenitor cells in the FVB/N-RC mice.
        PLoS One. 2012; 7: e43624
        • Díaz-Cruz E.S.
        • Sugimoto Y.
        • Gallicano G.I.
        • Brueggemeier R.W.
        • Furth P.A.
        Comparison of increased aromatase versus ERα in the generation of mammary hyperplasia and cancer.
        Cancer Res. 2011; 71: 5477-5487
        • Frech M.S.
        • Halama E.D.
        • Tilli M.T.
        • et al.
        Deregulated estrogen receptor alpha expression in mammary epithelial cells of transgenic mice results in the development of ductal carcinoma in situ.
        Cancer Res. 2005; 65: 681-685
        • Fan L.
        • Zhang F.
        • Fan H.
        • Zhang C.
        Brief review of image denoising techniques.
        Vis Comput Ind Biomed Art. 2019; 2: 7
        • Frech M.S.
        • Torre K.M.
        • Robinson G.W.
        • Furth P.A.
        Loss of cyclin D1 in concert with deregulated estrogen receptor alpha expression induces DNA damage response activation and interrupts mammary gland morphogenesis.
        Oncogene. 2008; 27: 3186-3193
        • Kumar N.
        • Nachamai M.
        Noise removal and filtering techniques used in medical images.
        Orient J Comput Sci Technol. 2017; 10: 103-113
        • Meijering E.
        • Jacob M.
        • Sarria J-Cf
        • Steiner P.
        • Hirling H.
        • Unser M.
        Design and validation of a tool for neurite tracing and analysis in fluorescence microscopy images.
        Cytometry A. 2004; 58: 167-176
        • Tolg C.
        • Cowman M.
        • Turley E.A.
        Mouse mammary gland whole mount preparation and analysis.
        Bio Protoc. 2018; 8: e2915
        • van der Walt S.
        • Schönberger J.L.
        • Nunez-Iglesias J.
        • Boulogne F.
        • Warner J.D.
        • Yager N.
        • Gouillart E.
        • Yu T.
        • The scikit-image Contributors
        scikit-image: image processing in Python.
        PeerJ. 2014; 2: e453
        • Harris C.R.
        • Millman K.J.
        • van der Walt S.J.
        • et al.
        Array programming with NumPy.
        Nature. 2020; 585: 357-362
        • Hunter J.D.
        Matplotlib: a 2D graphics environment.
        Comput Sci Eng. 2007; 9: 90-95
        • Virtanen P.
        • Gommers R.
        • Oliphant T.E.
        • et al.
        • SciPy 1.0 Contributors
        SciPy 1.0: fundamental algorithms for scientific computing in Python.
        Nat Methods. 2020; 17: 261-272
        • Sherratt M.J.
        • McConnell J.C.
        • Streuli C.H.
        Raised mammographic density: causative mechanisms and biological consequences.
        Breast Cancer Res. 2016; 18: 45
        • Nazari S.S.
        • Mukherjee P.
        An overview of mammographic density and its association with breast cancer.
        Breast Cancer. 2018; 25: 259-267
        • Mainprize J.G.
        • Alonzo-Proulx O.
        • Alshafeiy T.I.
        • Patrie J.T.
        • Harvey J.A.
        • Yaffe M.J.
        Prediction of cancer masking in screening mammography using density and textural features.
        Acad Radiol. 2019; 26: 608-619
        • Kerlikowske K.
        • Scott C.G.
        • Mahmoudzadeh A.P.
        • et al.
        Automated and clinical breast imaging reporting and data system density measures predict risk for screen-detected and interval cancers: a case-control study.
        Ann Intern Med. 2018; 168: 757-765