Astrocyte regional heterogeneity revealed through machine learning-based glial neuroanatomical assays

Abstract

Evaluation of reactive astrogliosis by neuroanatomical assays represents a common experimental outcome for neuroanatomists. The literature demonstrates several conflicting results as to the accuracy of such measures. We posited that the diverging results within the neuroanatomy literature were due to suboptimal analytical workflows in addition to astrocyte regional heterogeneity. We therefore generated an automated segmentation workflow to extract features of glial fibrillary acidic protein (GFAP) and aldehyde dehydrogenase family 1, member L1 (ALDH1L1) labeled astrocytes with and without neuroinflammation. We achieved this by capturing multiplexed immunofluorescent confocal images of mouse brains treated with either vehicle or lipopolysaccharide (LPS) followed by implementation of our workflows. Using classical image analysis techniques focused on pixel intensity only, we were unable to identify differences between vehicle-treated and LPS-treated animals. However, when utilizing machine learning-based algorithms, we were able to (1) accurately predict which objects were derived from GFAP or ALDH1L1-stained images indicating that GFAP and ALDH1L1 highlight distinct morphological aspects of astrocytes, (2) we could predict which neuroanatomical region the segmented GFAP or ALDH1L1 object had been derived from, indicating that morphological features of astrocytes change as a function of neuroanatomical location. (3) We discovered a statistically significant, albeit not highly accurate, prediction of which objects had come from LPS versus vehicle-treated animals, indicating that although features exist capable of distinguishing LPS-treated versus vehicle-treated GFAP and ALDH1L1-segmented objects, that significant overlap between morphologies exists. We further determined that for most classification scenarios, nonlinear models were required for improved treatment class designations. We propose that unbiased automated image analysis techniques coupled with well-validated machine learning tools represent highly useful models capable of providing insights into neuroanatomical assays.

Keywords: astrocyte; clustering analysis; gliosis; machine learning; neuroanatomy; neuroinflammation.

Figure 1.. Traditional Image Analysis Methods

Immunofluorescence imaging was used to determine labeling of astrocytes with GFAP and ALDH1L1. Images from vehicle and LPS-treated mice were used to determine the pixel intensity and colocalization of GFAP and ALDH1L1. A scatterplot was generated with a random sample of 20% of the data from each region (data shown only for hippocampus, a1/a2). There is no apparent difference in the colocalization of GFAP and ALDH1L1 pixel intensity in vehicle or LPS-treated mice (a3/a4). Similar plots were observed for neocortex and caudoputamen. Next, evaluation of the mean fluorescence intensity was calculated for GFAP (b1) and there was no significant increase between the vehicle and LPS (Hippocampus- p = 0.39; Neocortex- p = 0.47; Caudoputamen- p = 0.41). ALDH1L1 mean pixel intensity was calculated next (b2) and there was no difference between the control and treated groups (Hippocampus- p = 0.44; Neocortex- p = 0.43; Caudoputamen- p = 0.82).

Figure 2.. Image Capture Methodology

Images were captured at 40x magnification of three regions on the same coronal section of 5 vehicle-treated and 5 LPS-treated mice. Images in the caudoputamen and hippocampus were obtained ipsilaterally while the neocortex images were captured bilaterally on the same section. Representative images of each region are shown (a1/a3-b1/b3). 9-12 images were captured per mouse in each region with 152 images total per condition.

Figure 3.. Otsu Segmentation Workflow

Immunofluorescence imaging was used to determine labeling of astrocytes with GFAP and ALDH1L1. To evaluate accuracy (a1), specificity (a2), and sensitivity (a3), segmented masks were evaluated against two ground truth images. Otsu outperformed the other segmentation methods. The Otsu workflow is described in b and c. Images were average projected in FIJI and loaded into R studio. Separate the channels into GFAP positive and ALDH1L1 positive labeled astrocytes (b1/c1). Create a grayscale image with pixel intensities from 0-255 (b2/c2). Apply filter and Otsu threshold to reduce the gray-level image to a binary image with intensities 0 or 1 (b3/c3). Process image to fill holes between pixels of a specified distance. Apply segmented boundary to original image and ensure object was properly segmented. Apply gamma correction as needed to reduce background luminance levels. Compute features of each object to extract morphological and textural features for analysis between conditions and anatomical regions of interest (b4/c4). (d1)- 94.45% accuracy of random forest model predicting ALDH1L1 and GFAP segmented objects obtained from the hippocampus of vehicle-treated mice (p < 2.2 x e⁻¹⁶ over the NIR). (d2)- Random forest model created with GFAP vehicle-treated objects from all regions with overall accuracy of 79.26% and balanced accuracy of 94.49%, 78.95%, 75.72% for the hippocampus, neocortex, and caudoputamen, respectively. Significance is based on calculating the probability that the accuracy elicited from the model is higher than the NIR. (d3)- Random forest model created with ALDH1L1 vehicle-treated objects from all regions used to predict the region in which objects originate. Overall accuracy is 84.07%. Balanced accuracy is 93.32%, 83.84%, 85.73% for hippocampus, neocortex, and caudoputamen respectively. (e1-e3)- Average of categories for Gini Importance in random forest models displayed in d1-d3 respectively. Pixel intensity features had the highest importance in each model. Specific model details found in supplemental table 2.

Figure 4.. Feature Reduction and Random Forest Models Predicting Reactivity State

Using objects extracted from the hippocampus, a random forest was generated to perform the Boruta algorithm to identify key variables. The hippocampal random forest reactivity model (a) was applied to the full compute.features hippocampal test set, neocortex, and caudoputamen and evaluated on performance of accurate prediction. The GFAP hippocampal model had an accuracy of a 78.28% on hippocampal objects and a significance of p < 2.2 x e⁻¹⁶. Applied to the neocortex and caudoputamen, it had an accuracy of 49.80% and 52.65%, respectively (p = 1). The ALDH1L1 hippocampal model had an accuracy of 78.69% (p < 2.2 x e⁻¹⁶) on hippocampal test data. On neocortex and caudoputamen objects, the accuracy was 41.71% and 58.09% (p = 1 and p < 2.2 x e⁻¹⁶, respectively). Twelve variables were selected based on a correlation plot of all features and other studies performing image analysis with Haralick features. New random forest models created with hippocampal objects identified each feature with Boruta algorithm (b) for importance. Correlation plot of twelve variables (c) used in the model based on 10000 objects from the hippocampus with each biomarker. Dark blue indicates high positive correlation, dark red indicates high negative correlation. Model performance listed in supplemental table 3.

Figure 5.. Regionalization Predicting Activation and Eigen Vectors Showing Drivers of Variation of GFAP Segmented Objects

Random forests models were created with 12 variables to predict reactivity state of segmented objects from one region. Models were applied to other regions and evaluated on performance of accurate prediction of reactivity state. Random forest models were accurate at only predicting object reactivity state in their regions (a1-c1), indicating a heterogenous response to inflammation (see supplemental table 3 for model performance statistics). To identify the drivers of variance between vehicle segmented objects and LPS segmented objects, a PCA was generated and eigenvectors calculated. (a2/b2) Vehicle segmented objects from the hippocampus and neocortex placed higher importance on shape features (area and perimeter) and lower importance on pixel intensity. In contrast, LPS segmented objects placed higher importance on pixel intensity and lower importance on shape features (a3/b3). In both vehicle segmented and LPS segmented objects from the caudoputamen, the highest driver of variance was Haralick features.

Figure 6.. Regionalization Predicting Activation and Eigen Vectors Showing Drivers of Variation of ALDH1L1 Segmented Objects

Random forests models were created with 12 variables to predict reactivity state of segmented objects from one region. Models were applied to other regions and evaluated on performance of accurate prediction of reactivity state. Random forest models had higher accuracy at predicting object reactivity state in their specific regions (a1-c1), however, models from the hippocampus and caudoputamen were statistically significance at predicting one another’s activation state indicating a similar response to inflammation (see supplemental table 3 for model performance statistics). To identify the drivers of variance between vehicle segmented objects and LPS segmented objects, a PCA was generated and eigenvectors calculated. (a2/b2) Vehicle segmented objects from the hippocampus and neocortex placed higher importance on shape features (area and perimeter). In contrast, LPS segmented objects placed higher importance on pixel intensity (a3/b3). In vehicle segmented objects from the caudoputamen (c2), the highest driver of variance was pixel intensity while LPS segmented objects (c3) placed higher importance on shape.

Figure 7.. Clustering Methods using K-means, K-medioid, DBSCAN, OPTICS, and Distribution-based

In the GFAP LPS object segmented from the hippocampus, multiple clustering methods were performed in order to elucidate the different ways objects inherently cluster. Using the silhouette method, k-means determined the optimal number of groups was k = 5. In order to see if these groups clustered similar ways using t-SNE, plots were created without including the cluster designation in the modeling dataframe. The PCA (a1) is color coordinated for the cluster assignments and a t-SNE plot (a2) shows that these objects have very little overlap other than where each cluster meets another. K-medoid shows similar performance to k-means and clusters generally do not overlap except for where one cluster is plotted next to the other (b1/b2). Next, DBSCAN determined the optimal number of clusters was 7 (minimum points = 30). Black objects are ones deemed noise and not included in any cluster nor do they create their own cluster. Predominately, objects fell into cluster 1 shown in indigo (c1/c2). Utilizing OPTICS, the optimal number of clusters was also 7 (minimum points = 30). Compared to the DBSCAN method, OPTICS only classified 2 objects differently in which they were moved from cluster 7 to the noise category (d1/d2). Although these clusters may form inherently through DBSCAN and OPTICS, they provide very little insight into the data due to the diversity of 20,511 objects classified into cluster 1. Distribution-based clustering was performed and determined the optimal number of clusters was 8. When plotting color-coded for the cluster designation, we can appreciate there are undefined cluster borders (e1/e2).

Figure 8.. Clustering is Meaningful

To determine the optimal number of clusters using 30 different clustering techniques, we utilized the NbClust package in R. In the GFAP LPS hippocampus objects, the optimal number was 3 clusters (a), but when running multiple iterations, cluster size varied and objects with distinct morphological profiles were grouped into one cluster further limiting our analysis (b; compare cluster designation to Figure 7 a1). Using the silhouette method (c) we determined the optimal number of clusters and ran a permutation analysis (d) to determine if our clusters were distinct compared to randomly assigned clusters. After 1000 iterations, the lack of overlap in the average distance between or within clusters and the within cluster sum of squares suggests the probability of this occurring due to chance as p < 10⁻² (e).

Figure 9.. Evaluation of Hippocampal Astrogliosis Clusters in GFAP

Within the hippocampus, although there were five clusters predicted with the silhouette method, only Cluster 1 and Cluster 3 included identifiable processes. Cluster 1 LPS objects are colored indigo in the t-SNE plot (a). Gray line intersection indicates the location of the LPS object (c) shown in relation to other objects within the cluster in the t-SNE plot. Generally, these were large objects (~5 μm) and contained the cell body and long processes (b/c). Cluster 3 objects are colored in teal in the t-SNE plot (a). These objects were generally segmented possesses in cross section or a small portion longitudinally (~ 1 μm, b). When generating a LDA, there is a large overlap between the vehicle and LPS objects (d). The most important variables classifying vehicle versus LPS-treated objects were fluorescent intensity and texture features (e). Although there is large overlap between the objects as seen in the LDA density plot, the random forest from figure 5 (a1) was able to predict with high accuracy the reactivity of objects (Cluster 1- Region Specific (RS): 73.53%; Cluster Specific (CS): 72.75%. Cluster 3- RS: 66.62%; CS: 65.69%).

Figure 10.. Evaluation of Hippocampal Astrogliosis Clusters in ALDH1L1

In ALDH1L1 objects, only Clusters 1, 4, and 6 contained recognizable objects. Gray line intersection indicates the location of the LPS object (c) shown in relation to other objects within the cluster in the t-SNE plot (a). Examples of vehicle and LPS objects are shown (b/c). Similarly to the GFAP labeled objects, the LDA shows a high overlap of vehicle and LPS segmented objects but the random forest was highly accurate and every model had a significance of p < 2.2 x e⁻¹⁶. Accuracies are Cluster 1- Region Specific (RS): 78.25%; Cluster Specific (CS): 73.61%; Cluster 4- RS: 88.45%; CS: 79.91%; Cluster 6- RS: 74.46%; CS: 74.25%.