Introduction

In the healthy brain, microglial cells are characterized by a complex arborization of highly ramified and fine processes that are evenly distributed around the cell body. Their processes present a strong and constant dynamic of protractions and retractions. For many years, this dynamic has been associated with their monitoring role, as they are constantly looking for danger signals in their surroundings (Davalos et al., 2005; Nimmerjahn, Kirchhoff, & Helmchen, 2005). Microglia continuously adapt to their local environment, sensing even small changes via membrane receptors and transporters (Hanisch & Kettenmann, 2007). They are able to react quickly and adjust their functions in response to the detected signals. These functional changes are supported by morphological modifications. As an example, the formation of filopodia enables directed migration to lesion sites, while the formation of phagocytic cups and phagolysosome increases are linked to their phagocytic activity (Kettenmann, Hanisch, Noda, & Verkhratsky, 2011; Sierra, Abiega, Shahraz, & Neumann, 2013).

The traditional view of microglia being either ramified (homeostatic microglia) or amoeboid (reactive microglia) has now been proven to be oversimplified and inaccurate, as microglia display a whole range of various morphologies (Paolicelli et al., 2022). Indeed, many studies have shown a wide diversity of microglial morphotypes in pathological conditions, and microglia can even exhibit an increase in ramification within the first hours of injury (Vidal-Itriago et al., 2022).

In stroke models, morphological transformation of reactive microglia is correlated with the severity and duration of ischemia, as well as an enhanced expression of markers such as CD11B and CD68 (Zhang, 2019). In the rodent brain, different morphotypes can be observed depending on how the area has been impacted by the metabolic challenge. Thus, a gradient of microglial morphologies can be observed from amoeboid microglia in the ischemic core to ramified complex cells in the contralateral cortex (Anttila, Whitaker, Wires, Harvey, & Airavaara, 2017; Fumagalli, Perego, Ortolano, & De Simoni, 2013). Similarly, in Alzheimer’s disease (AD) mouse models, microglial cells associated with Aβ plaques present strong morphological transformations in opposition to plaque-distant cells mildly morphologically impacted (Plescher et al., 2018).

As a matter of fact, morphology is an indicator that is easy to obtain when characterizing a model through live imaging or immunohistochemistry, and it has been used in numerous studies to describe a broad range of pathological contexts (Adeluyi et al., 2019; Ali et al., 2019; Morrison, Young, Qureshi, Rowe, & Lifshitz, 2017; Plescher et al., 2018). To date, the literature contains a wide variety of criteria to quantitatively describe microglial morphology, ranging from descriptive measures such as cell body surface area, perimeter, and process length to indices calculating different parameters such as circularity, roundness, branching index, and clustering or artificial intelligence (AI) approaches to categorize morphologies (Adaikkan et al., 2019; Heindl et al., 2018; Kongsui, Beynon, Johnson, & Walker, 2014; Leyh et al., 2021; Morrison et al., 2017; Young & Morrison, 2018). This variety of metrics and approaches to measure microglia morphology makes it difficult to compare results between studies. Moreover, categorization methods seem subjective with the number, name and composition of categories being non consistent between studies. In addition, as microglial morphology is highly diverse, it is not rare to be confronted with cells belonging to categories that are not defined, or that could belong to several of them.

Here, we propose a fully automated program that quickly and unbiasedly computes most parameters used in the literature to characterize microglial morphology, and overcomes the limitations of current methods. Our approach consists in: 1) performing principal component analyses (PCA) to select the parameters best fitted to describe the dataset of interest and 2) generating a sorting of the cells based on their morphology. The originality of our study is the use of Andrews Plots to accurately rank microglia according to their morphology and discriminate distinct populations. In addition, our working flow has been developed so that it can be used for a wide variety of models and image acquisitions (confocal, wide field, two-photon, etc.). We developed this method on an exploratory dataset from a rat brain with focal cerebral ischemia, and validated it on other independent datasets from other models (other animal, age, pathology) to confirm its wide range of application.

Methods

Models

All experimental procedures were carried out in accordance with the French institutional guidelines and ethical committee, and authorized by the local Ethics Committees of the University of Bordeaux and “Comité d’éthique pour l’Expérimentation Animale Neurosciences Lyon”: CELYNE; CNREEA no. 42, and the Ministry of National Education and Research (APAFIS#30666-2021032518063894).

Ischemic stroke rat model – fixed tissue

The experiment was conducted on a 6-8-week-old Sprague Dawley rat with transient (90 minutes) middle cerebral artery occlusion (tMCAo) followed by reperfusion. Twenty-four hours after reperfusion, the animal was deeply anesthetized (xylazine at 12 mg/kg and ketamine at 90 mg/kg in 0.9% NaCl) and died by the collapse of the lungs when the thoracic cage was open. Infusion of phosphate-buffered saline (PBS) was performed directly in the left ventricle. Next, 100 mL of paraformaldehyde (PFA) 4% was infused at 10 mL/min before the brain collection. The brain was then post-fixed in PFA 4% for 2-4 days, rinsed, and soaked in 30% sucrose for 72 hours before being frozen in −40°C 2-methylbutane. It was then stored at −20°C until being cut. Floating 30 µm thick sections were obtained using an NX50 cryostat and conserved in PB-Thimerosal until immunostaining.

Ischemic stroke mouse model – live tissue

An ischemic stroke was induced in 8 to 12 weeks CX3CR1^+/GFP mice by a permanent occlusion of their middle cerebral artery (pMCAo). Imaging was performed the next day directly on living animals through a cranial window. Briefly, animals were anesthetized with isoflurane (induction: 3-4%, surgery: 1.5-2%) and installed on a stereotaxic apparatus where their temperature was monitored and maintained at 37°C. After cleaning and exposition of the skull, a polyamide implant was glued on a zone at the periphery of the lesion in order to be able to image both ischemic and extralesional regions. The skull was then cautiously thinned by drilling until a thickness of 20 to 30 µm was reached. A glass coverslip was then glued to the thinned region. Mice were then put under a heating lamp and individually housed until imaging (V. Hubert et al., 2021).

Phox2b mutation and Embryonic dataset

The expansion of a 20-residue polyA tract in PHOX2B is used as a model for Central Hypoventilation Syndrome (CCHS) (Amiel et al., 2003). Mutant pups were generated by mating conditional Phox2b^27Ala/27Ala females with Pgk::Cre males (Dubreuil et al., 2008). Since mutant newborns die rapidly after birth, all experiments were performed blind on E18.5 embryos of either sex, and the genotype of each embryo was determined a posteriori on tail DNA. Pregnant mice were killed by cervical dislocation at E18.5. Embryos were rapidly excised from uterine horns and placed in artificial oxygenated cerebrospinal fluid (aCSF) at 18-20°C until dissection. The aCSF solution composition was (in mM): 120 NaCl, 8 KCl, 0.58 NaH2PO4, 1.15 MgCl2, 1.26 CaCl2, 21 NaHCO3, 30 Glucose, pH 7.4 and equilibrated with 95% O₂-5% CO₂. Embryos were decerebrated and decapitated, and ventral tissues were removed. Brainstems were then carefully disassociated from the surrounding tissues and completely isolated by a rostral section made at the junction between the rhombencephalon and mesencephalon and a caudal section at the spinal cord level. After dissection, brainstem preparations were placed in 4% PFA for 2 - 3 h for tissue fixation. To obtain transverse frozen sections, brainstems were placed in a 20% sucrose-PBS solution for cryoprotection overnight, then they were embedded in a block of Tissue Tek (Leica Microsystems, France) and sectioned at 30 µm using a cryostat (Leica Microsystems, France).

Alzheimer’s Disease model

C57BL6 5-month-old female APPxPS1-KI and control (PS1-KI) mice were used²⁵. For the brain collection, 0.5 mg xylazine was injected intraperitoneally followed by Euthasol. They were then intracardially perfused first with 4% PFA in PB. After one hour, the brains were collected and then post-fixed for 2 hours in 4% PFA and stored in PBS. They were rinsed in phosphate buffer (PB) and soaked in 30% sucrose for 48 hours before being frozen at −40°C in 2-methylbutane for 2 minutes and stored at −70°C until the cut. Free-floating sections (30 µm thick) were cut using a cryostat and conserved in PB-Thimerosal until immunostaining.

In vitro microglia

We followed the protocol published by Bohen and collaborators using Magnetic Activated Cell Sorting (MACS) (Bohlen, Bennett, & Bennett, 2019). Ten days old mice were anesthetized with isoflurane, intracardially perfused with PBS enriched with magnesium and calcium (14040141, Gibco) before decapitation and brain collection. The brains were then finely chopped in a dissociation buffer and transferred to a tissue grinder to dissociate further the cells from the tissue. Myelin and debris were then eliminated thanks to a Percoll® PLUS solution (E0414, Sigma-Aldrich) and an incubation with anti-myelin magnetic beads (130-096-433, Miltenyi Biotec): the cell suspension is applied on a column fixed to a magnet so that the myelin will stay on the column while the cells will exit. To isolate specifically microglia, the suspension was incubated with anti-CD11b magnetic beads (130-097-142, Miltenyi Biotec). Cells were then seeded on 35 mm diameter dishes with 4 compartments with a polymer bottom (80416, Ibidi) coated with poly-D-lysine (100 µg/mL, A-003-E, Sigma-Aldrich) and collagen IV (2µg/mL, 354233, Corning). The culture medium was adapted from Bohlen et al., 2019 with DMEM/F12 (21041025, Gibco), 1% penicillin-streptomycin, 2 mM L-Glutamine (10378016, Gibco), 5 μg/mL N-acetylcysteine (A9165, Sigma-Aldrich), 100 μg/mL apo-transferrin (T1147, Sigma-Aldrich), 100 ng/mL de sodium selenite (S5261, Sigma-Aldrich), 1.5 μg/mL cholesterol (700000P, Avanti), 0.1 μg/mL oleic acid (90260, Cayman Chemical), 1 ng/mL gondoic acid (20606, Cayman Chemical), 1 μg/ml heparan sulfate (AMS.GAG-HS01, AMSBIO), 2 ng/mL TGFβ2 (100-35B, Preprotech) and 10 ng/mL M-CSF (315-02, Preprotech). Cells were kept at 37°C 5% CO₂ and the culture medium was half renewed every two days until imaging (4 to 10 days after seeding).

Immunohistochemistry

The sections underwent 5-minute rinses three times with PBS and Triton 100X 0.3% (PBS-T). Sections from the ischemic experiment were incubated for 20 minutes in a 1.5% hydrogen peroxide solution diluted in PBS, and after 3*10 minutes of rinsing with PBS-T, they were incubated for 30 minutes in 0.5% sodium borohydride diluted in PBS and then thoroughly rinsed. One-hour saturation was then performed on all sections with PBS-T, 3% bovine serum albumin and 5% normal goat serum (PBS-T-BSA-NGS). Sections were then incubated overnight with anti-Iba1 antibodies (1/500, Wako, Germany, #W1W019-19741) at 4°C and anti-Amyloid-β 4G8 (1/500, Bio-Legend, #800709) antibody. Sections were rinsed 3*10 minutes the next day with PBS-T. They were then incubated in goat anti-rabbit secondary antibodies coupled with Alexa 555 fluorochrome (1/300, Invitrogen, France, #A21429) and goat anti-mouse coupled with Alexa 488 fluorochrome (1/300, Invitrogen, France, #A11029) for 2 hours at room temperature. They were then rinsed 3*10 minutes with PBS and mounted on microscope slides using Fluoromount G (ThermoFisher Scientific, #00-4958-02) and a glass coverslip.

For the Phox2b mutant model and embryonic dataset, sections were incubated for 90 min in a solution of PBS containing 0.3% Triton X-100 and 1% BSA to limit nonspecific labeling and favor antibody tissue penetration. The primary rabbit anti-Iba1 antibody (1/500; Wako, Germany, #W1W019-19741) was then applied overnight at room temperature and under slight agitation. To amplify the primary antibody signal, after several washes, sections were incubated with an Alexa Fluor 568 donkey anti-rabbit secondary antibody (1/500, Merck, France, #A10042) for 90 min at room temperature. Immunostained sections were then mounted in Vectashield Hard Set medium (Eurobio, France) cover-slipped and kept in the dark until imaging.

Image acquisitions

For the ischemic stroke model on fixed tissue, acquisitions were performed with a Zeiss confocal microscope and a 40X objective (LSM 880, Zeiss). Z-stacks between 18 and 25 µm thickness were executed with a 1 µm step. For the ischemic stroke model on live, microglia were observed thanks to a two-photon microscope (Bruker Ultima) and a 20X objective, at a deepness of 50 to 150 µm. 10 to 15 µm thick acquisitions were acquired with a 1 µm step between plans. For the AD model, a Zeiss Axio Scan 7 microscope slide scanner was used with a 20X objective. Z-stacks (30 µm thick) were performed with a step of 1 µm. For the Phox2b mutant model, imaging was conducted using a Zeiss confocal microscope (LSM 900, Zeiss). Z-stacks (30 µm thick) were performed with a step of 1 µm and a 40X objective. In vitro microglia were imaged on a single plane using the BioStation (Nikon).

Image processing and segmentation

The general procedure is described in supplemental Fig. 7.

For the ischemic model on fixed tissue, all post-acquisition processing was executed with FIJI software. Individual microglial cells were cropped from the original acquisitions before Z-maximum projections. The brightness of the projections was adjusted, and median (radius=2), mean (radius=2) and unsharp mask (radius=5; mask=0.60) filters were applied via a macro. A threshold was then applied to each crop to binarize the cells, and pixels not belonging to the microglia of interest were manually deleted. The same procedure was performed for the Phox2b mutant model, except the mean filter that was not applied.

For the ischemic model on live tissue all post-acquisition processing was also executed with FIJI software. Z-maximum projections were performed, and individual crops of the cells were done. Brightness and contrast of the images were modified before applying a manual threshold to binarize the cells. Reminiscent noise and neighboring cells were manually removed.

For the AD model, crops of individual microglial cells located in the secondary visual cortex were generated with Zen software and exported to the Tif image format. The next operations were executed with FIJI software. Z stack cleaning was performed by background subtraction (rolling=50), followed by a median filter (radius=3) and an unsharp mask (radius=5; mask=0.60). Segmentation was performed by the FIJI Trainable Weka Segmentation 3D plugin (Arganda-Carreras et al., 2017). The selected training features were based on the Gaussian blur, Laplacian, Maximum, Mean, Minimum, Median and Variance (Minimum sigma: 1.0 and Maximum sigma: 8.0). Binarized z-stacks were then manually cleaned to remove pixels from neighboring cells or remaining noise. Finally, a maximum Z-projection was performed.

For the Phox2b mutant model, image processing was performed via a custom MATLAB R2018B (Mathworks, Massachusetts, USA) script. First, a median filter followed by an unsharp mask filter was applied to all images. An automated threshold was then applied. Manual tuning with a thresholding tool (Bemis, 2023) was applied if the automated threshold was not satisfying. The binarized images were then automatically cropped around microglia to manually remove pixels belonging to neighboring cells.

For the in vitro microglia, acquisitions were cropped around individual microglia, and the cells were then manually traced on Procreate (5.2.9 version).

When the binarization process generated one or several holes inside the cell bodies, they were manually filled to not interfere with the automated detection of the cell bodies. All binarized and individualized microglial cells were then rescaled through a FIJI macro command to obtain a pixel resolution equal to 0.5 µm. Finally, the images were resized to obtain a field size equal to 300 x 300 pixels.

MorphoCellSorter algorithm pipeline

Automated measure of 20 dimensionless morphological descriptors

All parameters are illustrated in Supp. Fig. 8.

First, we calculate the morphological parameters most commonly used in the literature, followed by other parameters that we have developed and proposed here (Fig. 2B). To make our approach scale-independent, we reformulated some morphological descriptors to obtain dimensionless parameters.

Cell Body Recognition

The cell body identification is based on successive erosions until the center of the cell body is the only pixels remaining using a 3*3 pixels square as a structuring element. The last erosion is the one preceding the elimination of all pixels belonging to the cell and corresponds to the center of the cell body. Then, successive dilations are performed to reconstruct the cell body according to the morphology of the cells, using a cross (+) 3 pixels height and width (Leyh et al., 2021). Reconstruction of the cell body relies on (1) the identification of the index i on which measured surfaces after erosions are inferior to the surfaces measured after dilation, (2) the average of two copies of the images: one on which i + 1 erosions have been performed, and the other one on which i dilatations have been performed from the center of the cell body. The resulting mean image was then subjected to a threshold fixed at 0.25 to obtain the cell body. From this, we can calculate the following parameters.

Skeleton

The cell skeleton is obtained via the MATLAB bwskel function. This function reduces all 2-D binary objects to 1-pixel width curves while taking care of not changing the essential structure of the image. The aim is to extract the image skeleton while keeping the topology properties such as the Euler number and to extract branchpoints and endpoints. The skeleton length is defined by the number of pixels forming the skeleton minus the pixels of the cell body.

Thus, we calculate the skeleton related parameters:

Sholl analysis

Sholl analysis was performed on the microglial skeleton. Concentric circles with a step size of two pixels were generated on the skeletons minus the soma, and the soma centroid was the center point of the circles. Intersections between the growing circles and the skeleton were counted.

Fractal/Hausdorff dimension

The fractal dimension measures the complexity of a cell shape, that is, how an object fills the space at different scales. The Hausdorff fractal dimension (Costa, 2023) corresponds to the slope of the logarithm of the number of tiles of size ε needed to cover all the surfaces of the object of interest vs the logarithm of ε.

Lacunarity

Lacunarity corresponds to the inhomogeneity of a cell shape. It measures the variance of how the object fills the space at different scales. Patterns with larger gaps generally have higher lacunarity. We used the gliding box algorithm (Vadakkan, 2023; Tolle, McJunkin, & Gorsich, 2008), and the value corresponds to the slope of the log-log representation.

Convex Hull

The convex hull is the smallest convex polygon (with angles inferior to 180°) enclosing the whole cell.

We also measure the length of the major and minor axes of the convex hull.

Finally, we measure the largest and smallest radii of the circles from the centroid of the convex hull as the center to external points.

Linearity

A Principal Component Analysis is generated on all the points forming the microglial cell. Then, the variance of the two first principal components (PCs) is calculated; the highest corresponds to the Variance Max and the lowest to the Variance Min.

Establishment of the cell ranking according to morphological criteria

The second step is to determine which of the previously calculated parameters best discriminate the cells in a given dataset. Principal component analysis (PCA) is a statistical method allowing dataset dimensionality reduction (Jolliffe & Cadima, 2016). A PCA is performed on the 20 morphological indexes calculated to determine the parameters that best explain the dispersion of the data. Lacunarity, roundness factor, convex hull radii ratio, processes cell areas ratio and skeleton processes ratio were subjected to an inversion operation in order to homogenize the parameters before conducting the PCA: small values for less ramified or more linear cells and higher values for more complex cells less linear. To identify the main parameters responsible for the data spreading, we first transform the raw dataset table into a z-scored table by the transformation (1), where < X̃_j > and σ_x̃j are respectively the average and the standard deviation of the parameter F_j. Because the amplitudes of the parameters are generally very different, such a transformation is necessary to compare their ranges in the same referential.

We next calculate the parameters correlation matrix, as follow:

C_XX is a N × N symmetric square matrix (ρ_ij = ρ_ji, and ρ_ii = 1) which expresses the correlation between two parameters F_i, F_j taken in pairs. Note that | ρ_ij |≤ 1. This matrix is comparable to an inertia tensor of the data cloud. Thus, to find the main directions of the N dimensional data cloud, we calculate the eigenvalues and eigenvectors of C_XX, such that:

Solving this equation, we get then N eigenvalues λ₁,…,λ_N, each associated to their corresponding eigenvector V₁,…,V_N. Each of these eigenvectors has N components. The eigenvector associated with the biggest eigenvalue is the first principal vector or principal component PC₁, while the second principal component PC₂ is the eigenvector corresponding to the immediately smaller eigenvalue λ₂ < λ₁, and so on. Fig. 2B shows an example cascade of eigenvalues in descending order.

PC₁ is the direction in which data are most dispersed, that is the direction holding the most important part of the data information. PC₂ is the second direction holding the next important information part, and so on. The two first components hold almost the whole of information, hence we can consider the plan PC₁, PC₂ as the principal plan reducing the dataset to a two dimensional space. Next, we calculate the projection of the parameters F_n onto the two first principal components respectively PC₁ and PC2 (Fig. 2C) and rank those parameters according to the intensity (absolute value) of their projection onto PC₁ (Fig. 2D). Each parameter is then weighted by a factor W_n= w_n./w₀, where w_i is the absolute value of the n^th parameter projection onto PC₁, and w₀ is the absolute value of the strongest projection among all the parameters. As shown in Fig. 2C, each parameter F_n is represented by a vector in the PC₁, PC₂. Thus, we calculate the projection y_n of each parameter onto PC₂ and calculate the following phase factor ϕ_i = arctan(y_n/w_n).

The following stage of our method consists in introducing the previously calculated factors (Weight: W_n, and phase: ϕ_n) and their corresponding parameter into a special Fourier synthesis called “Andrews Plot” or “Andrews curve” (Andrews, 1972; García-Osorio & Fyfe, 2005) allowing to visualize the data structure in high-dimensional. The general form of the spectral Fourier synthesis of an individual j writes:

that we rewrite as,

where n, is the index of the ranked parameters, N is the full number of parameters F_n, t is a virtual time 0 ≤ t ≤ 1, and C_n = F_n × W_n. Finally, we derive this method to represent data in high dimensional by a simple one-dimension signal composed of orthogonal trigonometric functions in researching the time t = t^* for which the curves present the maximum of variance. At this specific time point, we rank individuals according to the magnitude of their own function S_j(t^*) value as represented Fig. 2E. We now call S(t^*) the “Andrews-Score”.

Statistics

Rankings comparisons, performed with Spearman’s correlation, were also used to measure the linear correlation between the rankings. They were executed with MATLAB.

We used mixed effect models on RStudio, that are adapted to complex experimental designs with dependent and non-dependent data. This allows the addition of random factors to fixed effects. Model fittings and analysis were performed with the lme4 package (Amiel et al., 2003). Condition (Control/APP-PS1xKI or WT/CCHS) was considered a fixed effect, and the factor mice was considered a random intercept. When the mixed models were not applicable, we used Wilcoxon tests. Statistical significance is shown on the graphs (*p < 0.05; **p < 0.01; ***p < 0.001).

Number of bins of Andrews scores histograms was automatically determined by the Freedman-Diaconis method based on the following formula:

With x_max and x_min being respectively the maximum and minimum Andrews scores, IQR the interquartile range and N the number of cells.

Results

Evaluation of classification methods

Morphological categorization has been at the center of recently developed approaches to characterize microglial morphology. It enables the evaluation of the morphological diversity existing in a cell population. Moreover, in the context of morphological comparison between populations, the study of the morphotypes’ proportions is relevant as it can highlight meaningful differences that the sole global comparison of morphological criteria cannot. One way to categorize the cells based on their morphology is to perform clustering to group the cells having similar characteristics. We tested a commonly used approach, k-means, to categorize microglia from a rat brain cortex following an occlusion of the middle cerebral artery (Fig. 1A). In this model, microglia contralateral to the lesion are mildly to not impaired whereas ipsilateral microglia undergo light to drastic morphological changes (Anttila et al., 2017; Fumagalli et al., 2013) (Fig. 1A). We selected 20 non dimensional parameters from the literature and developed by the team, on which we conducted a Principal Component Analysis (PCA). The k-means clustering algorithm was applied on the first two principal components (PCs), and we fixed at 4 the number of clusters as commonly chosen in the literature (Verdonk et al., 2016) (Fig. 1B).

Figure 1:

The obtained microglial groups could be named as ramified microglia (cluster 1), with cells having a morphology consistent with a physiological state. Cluster 2 could be referred to as “activated microglia”, as it was composed of microglia displaying various morphologies; they had an intermediary ramification level with large cell bodies. Cluster 3 was composed of cells having few short or no processes and corresponded to amoeboid microglia. Lastly, cluster 4 was rod-like microglia, with two processes extending on either side of the cell body or a single long polarized extension. All the cells of the dataset (ipsilateral and contralateral) and their associated cluster are represented in Fig. 1C. These 4 clusters were not equally represented across each condition as no amoeboid nor rod-like microglia were found in the contralateral side (Fig. 1D). The large majority of the cells were ramified microglia (92%), and the rest were activated microglia (8%). The proportion of ramified cells drastically decreased ipsilaterally by a reduction factor of 4 (23%), while activated microglia increased (39%). 25% of the cells were amoeboid microglia, while rod-like cells represented 13% of the ipsilateral population (Fig. 1D). These results show that the microglial populations in each region are different, with microglia having a non-pathological morphology mainly found contralaterally, with only few cells being potentially impacted by the ischemia/reperfusion challenge as they showed morphological signs of reactivity with a less ramified morphology. Ipsilateral microglia were morphologically highly altered, with a decrease in ramified cells, replaced by activated and amoeboid microglia, and rod-like microglia.

However, the clusters themselves were highly morphologically heterogeneous, making these quantifications problematic (Fig. 1C). For instance, Cluster 1 (ramified microglia), was composed of cells with branched processes but also included cells that were visibly affected by the ischemic lesion, exhibiting sometimes short, few, or less branched processes, as well as elongated cell bodies. These microglia might have been more appropriate in Cluster 2 (activated microglia). The composition of the latter is also questionable, with the presence of branched processes that would fit better in the first cluster, or cells with few very short extensions that resemble some cells found in Cluster 3. Several non-bipolar or non-unipolar cells are also found in Cluster 2. Some of these cells not belonging to their affiliated clusters were highlighted by a white cross on Fig. 1C. We tested to decrease or increase the number of clusters with no improvement (data not shown).

Other classification methods exist, relying on supervised artificial intelligence (AI) algorithms to identify specific morphotypes. These approaches can have high rates of success (96% for the method developed by Leyh and collaborators (Leyh et al., 2021)), meaning that they are in theory efficiently classifying cells. However, the problems of cells which do not fit in a defined category, or to several categories can be raised. For example, in our dataset we found that about 60% of cells would not belong to any category or belong to several categories. In this context, the way the networks are trained might be very subjective to the experimenter and could influence the results. Moreover, depending on the papers, semantic problems are also met, with similar morphologies having different names or on the contrary same names having different definitions (Choi et al., 2022; Leyh et al., 2021), making the comparisons between studies complicated.

We know now that microglial morphology is rather a continuum than closed categories (Paolicelli et al., 2022). That is why we propose a continuous ranking rather than a classification of the cells.

Development of a ranking approach based on Andrews plots to order rat fixed microglia in a context of ischemic stroke lesion

We developed MorphoCellSorter, a novel ranking method, on the same dataset of rat microglia in the context of an ischemic stroke lesion, on fixed tissue acquired via confocal imaging. This model is known to induce major morphological modifications in microglia (Violaine Hubert et al., 2021), enabling us to constitute a heterogeneous dataset, with varied morphologies ranging from amoeboid to hyper-ramified.

The method is summed up on Fig. 2A, and takes as inputs binarized and individualized cells form projected Z stacks. To generate the ranking, we first gathered 20 morphological indices, that we also call parameters, either from the literature (Clarke, Crombag, & Hall, 2021; Fernández-Arjona, Grondona, Granados-Durán, Fernández-Llebrez, & López-Ávalos, 2017; Madry et al., 2018) or developed by the team (see Methods). These parameters, describing morphological complexity (circularity, ramification index, roundness factor, perimeter area ration, density, processes soma area ratio, processes cell area ration, skeleton processes ratio, branchpoints endpoints ratio, fractal dimension, lacunarity, solidity, convexity, convex hull circularity, convex hull radii ratio and branching index), linearity (linearity, convex hull span ratio, inertia) and directionality (polarization index) of the cells, were normalized to suppress any dimension in order to be exempted from any scale effect. The first step of our method involves PCA, which helps us determine the main parameters responsible for the data dispersion. The first principal component (PC1) captures the direction of maximum variance in the data, and subsequent component (PC2) captures the direction of the next highest variance (Fig. 2B). Typically, the first two principal components account for a significant proportion (often over 70%) of the data’s total variance, allowing us to reduce the dataset’s complexity to a two-dimensional space, such as a two-dimensional plane defined by PC1 and PC2 (Fig. 2C).

Figure 2:

Next, each parameter was projected onto these principal components, which helps us rank parameters based on their influence on the overall data structure (Fig. 2D). To determine the number of parameters used to elaborate the ranking, we selected the median number of the normalized cumulative projection (threshold = 0.5). We assign weights to the selected parameters based on their projections onto PC1, which allows us to emphasize the most influential parameters. This weighted projection, along with phase factors calculated for each parameter (referred as ϕ_n in the Methods section), forms the basis of an Andrews Plot — a graphical representation that simplifies high-dimensional data into a two-dimensional signal using trigonometric functions.

The Andrews Plot represents each individual microglia, as a unique curve, allowing for visual comparison based on their trajectories within this synthesized signal (Fig. 2E). This method enables efficient visual comparisons and insights into complex datasets with potentially significant implications for data analysis and interpretation. We took advantage of this representation to rank the cells based on their morphology; the ranking is established at the time point displaying the highest variance, i.e. the time point where the curves are the most dispersed, also named 1 (Fig. 2E).

The obtained ranking is displayed on Fig. 3A. To evaluate the quality of the ranking method, we compared the automated ranking to a subjective manual ranking by experts based on visual morphological criteria (cell size; number, length and branching of the processes). To assess the similarity degree between the rankings, we calculated Spearman’s rank correlation coefficient (r_S=0.96; Fig. 3B,E). This test revealed that the automated ranking is close to the manual ranking. To ensure the reliability of the manual ranking, we compared the automated ranking with another expert and obtained the same results (r_S=0.97; Fig. 3C,F). Interestingly, the comparison between the two manual rankings shows similar results to those obtained when compared to the automated ranking (r_S=0.93; Fig. 3D,G), indicating that variability is the same between automated or manual rankings. However, automated ranking provided an objective and repeatable result, validating our method. Thus, we are able to generate an accurate ranking of the cells based on morphological criteria.

Figure 3:

Evaluation of MorphoCellSorter on other datasets of microglia from different animals, age, models and imaging techniques

One limitation to the use of a general and common method to characterize microglial morphology is the high diversity of study material: the models (in vivo, ex vivo, in vitro, fixed or live tissue), the type of collected signal (DAB or fluorescent staining, endogenous fluorescence…), the imaging technique (influencing the resolution and the background noise). We thus wanted to evaluate the capacity of MorphoCellSorter to accurately rank microglia from various datasets having different morphological characteristics (Table 1). As previously, we compared the results to manual rankings performed by two independent experts to evaluate the quality of the automated rankings (Table 1). We tested the accuracy of the rankings using several thresholds impacting the number of parameters selected as entries for the Andrews plots (Table 1). To maximize the accuracy of the automated ranking for all datasets, we determined which threshold provided the best rankings overall (Supp. Table 1). We observed that a threshold of 0.8 (Supp. Table 1), leading to 11 to 12 selected parameters (Supp. Fig. 1A-E), constituted a good compromise to obtain good rankings for all datasets. Thus, we chose to fix this threshold for the following analyses. The Andrews plots and rankings obtained at the 0.8 threshold are respectively displayed in Supplemental Fig. 2 and 3.

Table 1:

Thus, our Andrews plot-based ranking method, relying on an automatic selection of relevant parameters thanks to the PCA, allows a reliable automated ranking of the cells based on their morphologies. The diversity of the tested datasets, with different models, types of signal and imaging techniques shows that our approach offers a good sorting of the cells independently of their shape and context.

APPxPS1-KI microglia in the visual cortex have significant morphological alterations compared to control cells

To determine the usability of our approach to compare microglial morphology in different conditions, we used MorphoCellSorter to assess the presence of morphological alteration in the context of an adult pathological murine model of Alzheimer’s disease. This model is known to induce amyloid plaques with dystrophic neurites, synaptic dysfunction and behavioral deficits (Casas et al., 2004). We compared the morphology of PS1-KI (control) and APPxPS1-KI microglia in the visual cortex (Fig. 4A). Eleven parameters were automatically selected by the PCA: circularity, ramification index, perimeter area ratio, fractal dimension, skeleton processes ratio, density, convexity, processes cell areas ratio, lacunarity, branching index and processes soma areas ratio (Supp. Fig. 4A-C) and used to generate the Andrews plots (Supp. Fig. 4D). The automated ranking obtained revealed that control and APPxPS1-KI microglia were well segregated with APPxPS1-KI microglia at the beginning of the ranking and control mice at the end (Fig. 4B). The distributions showed a 16.4% overlap between the populations, meaning that a large number of APPxPS1-KI cells have a morphology that differs from that of controls (Fig. 4C).

Figure 4:

The morphological indices showed significant differences between APPxPS1-KI and control microglia. The parameters measuring the morphological complexity of the cells highlighted a clear decrease in the ramification level of APPxPS1-KI microglia, with for example a fractal dimension and convexity decrease and a circularity increase. However, there were no significant differences in the cells’ linearity and polarity between groups (Fig. 4D-W).

Thus, this approach allowed the detection of morphological changes of microglia in the visual cortex of APPxPS1-KI mice.

MorphoCellSorter pinpoints no morphological alterations in a model of congenital central hypoventilation syndrome

Congenital central congenital hypoventilation syndrome (CCHS) is a rare genetic disease is associated with mutations of the PHOX2B gene and characterized by life-threatening breathing deficiencies (Ceccherini, Kurek, & Weese-Mayer, 2022; Dubreuil et al., 2008). Phox2b is a transcription factor required for the specification of the autonomous nervous system that contains respiratory control centers. Because some of our preliminary observations strongly suggest a possible involvement of microglia in the disease (unpublished data), we tested whether MorphoCellSorter highlighted eventual morphological alterations of microglia in Phox2b mutants (CCHS embryos) compared to wildtype embryos (WT), focusing on brainstem tissue that hosts respiratory networks involved in breathing rhythmogenesis.

In this case using a threshold of 0.8, the PCA highlighted 12 parameters that were used to rank microglia with Andrews plots: circularity, ramification index, perimeter area ratio, skeleton processes ratio, fractal dimension, processes cell areas ratio, density, branching index, processes cell areas ratio, convexity, lacunarity and solidity (Supp. Fig. 5A-C). These parameters were weighted to generate the Andrews plots (Supp. Fig. 5D) used to rank the cells according to their morphology (Fig. 5B; complete ranking in Supp. Fig. 6).

Figure 5:

Both populations displayed a large variability of morphologies from amoeboid to few ramified processes and were well represented throughout the ranking (Fig. 5B). This is concordant with previous studies on developing microglia: they can concomitantly be found in a large variety of shapes and morphologies (Orłowski, Sołtys, & Janeczko, 2003). Indeed, the distributions of Andrews scores showed no difference between WT and mutant mice. The overlap of distributions between the populations found was 78.9% (Fig. 5C). Accordingly, when taking into consideration the various morphological indices, there was no significant difference between WT and Phox2b mutant microglial morphologies (Fig. 5D-W). Additional experiments are required to test whether this is specific for embryonic stages and whether this is associated or not to different microglial functional states.

Discussion

MorphoCellSorter enables the automated measurement of morphological descriptors, the calculation of several non-dimensional indices, and a meaningful comparison of cell populations based on their morphological characteristics. The method was initially applied to an ischemic stroke model on confocal images with significant morphological changes. MorphoCellSorter was subsequently tested on five different datasets with heterogeneous characteristics (various models, age, imaging technique and type of collected signal), showcasing its adaptability. We then successfully used our approach to compare the cells’ morphology in an Alzheimer’s disease context, and in embryonic microglia using a CCHS model. Our method can autonomously identify the most effective parameters for discriminating the studied cell set, thereby providing valuable insights into microglial morphology. Additionally, it offers a useful automated ordering of cells based on a combination of selected parameters using Andrews plots. This feature allows for precise positioning of cells along an axis, ranging from round cells to highly branched cells, facilitating the discrimination of various cell populations based on their morphologies.

The use of the Andrews plots offers the possibility to combine all automatically selected most discriminant parameters and to weight them according to their relevance, making the rankings more accurate and less subjective than only sorting the cells based on one parameter. Even if we recommend a threshold at 0.8 for the number of parameters selected, as it gave good rankings for all datasets we tested, we included the possibility to change this threshold if the generated ranking is not satisfying, knowing that all thresholds give accurate rankings when compared to expert rankings.

The automated parameter selection relies on principal component analysis, which is used to determine the parameters that best disperse the data. In contrast, other approaches use PCA for dimensionality reduction to create new morphological parameters that are difficult to apprehend as they do not refer to any physical measureclarke (Clarke et al., 2021; Heindl et al., 2018). The strength of MorphoCellSorter lies in its personalized approach, with automated parameter selection adapted to the dataset, making it easily applicable to segmented microglia in diverse contexts, including those with distinct morphologies, such as microglia in culture or in vivo acquisitions. It does not require program adaptation or the need to train a neural network for deep/machine learning approaches. We observed that some selected parameters were consistent between all the tested datasets (ramification index and perimeter area ratio for example), showing their robustness to differentiate microglial cells with various morphologies. On the contrary, some parameters were never selected (inertia and linearity for example). The information contained in these parameters were not discriminant enough in the datasets we used in this present work, but we can imagine that they would have been selected if we would have chosen a model in which microglial growth towards a specific site is observed, such as focal laser injury(Davalos et al., 2005; Haynes et al., 2006).

We also introduce a ranking tool to further enhance morphological analysis. In contrast to recent studies, we opt for ordering cells rather than classifying them. Microglial morphologies can span a continuum with no clear boundaries between morphotypes. Automated clustering methods such as k-means (Young & Morrison, 2018) may be more suggestive, as they require an a priori determination of the number of clusters. Hierarchical clustering (Clarke et al., 2021) can avoid this requirement but may result in non-reproducible clusters, or clusters without biological meaning. Attempts to identify distinct morphotypes, as previously done with artificial intelligence (Leyh et al., 2021), are also questionable due to varying definitions of morphotypes between laboratories and a lack of consideration for subtle morphological differences and intermediate forms. This problem is overcome by the ranking approach, as it does not lock the cells in defined categories but rather position the cells in relation to one another. However, the ability to compare the evolution of clusters of microglial morphologies in different conditions remains valuable, especially when morphologies are heterogeneous in the same population, such as in developmental microglia or pathological contexts. Comparing medians or means of certain parameters may not reveal differences due to high variability, whereas comparing the evolution of the proportions of specific morphotypes may provide more informative insights. The Andrews plots used to generate a ranking of the cells could also be used to generate clusters of the cells presenting morphological resemblances, as similar cells display similar Andrews curves’ trajectories.

We chose not to integrate a segmentation tool into our program, as the segmentation method should be tailored to the specific imaging conditions, materials, and microscopes used. Diverse segmentation methods are available in different software, with varying levels of automation and manual input requirements. In this work, we present methods using FIJI. It is important to note that achieving perfect cell segmentation is not a prerequisite for using MorphoCellSorter since the tool does not aim to provide highly precise measurements of microglial morphology. For instance, discontinuous processes are not problematic. Consequently, it is possible to work with lower-resolution microscopes, such as epifluorescence wide-field microscopes, although the segmentation process may be more labor-intensive if the source images are of lower quality. Obviously, the results of the ranking will be highly affected by the quality and the homogeneity of segmentation in a given data set. However, we did not assess how much the resolution of the microscope used affects the quality of the ranking as low resolution and poor signal to noise ratio would mainly affect the segmentation step which is not discussed in this paper.

One potential limitation is that we decided to work with 2D Z-stack projections, as accurate 3D reconstructions are time-consuming and require high-resolution acquisitions. Our goal was to develop a fast and robust method that could be broadly applicable to most microscopes and imaging scenarios. Therefore, our approach may not yield precise measurements of parameters such as the exact number and length of processes or their hierarchy. However, other research groups have successfully developed tools capable of providing these types of measurements, such as 3DMorph (York, LeDue, Bernier, & MacVicar, 2018) and MorphOMICS via Imaris (Colombo et al., 2022). Instead, MorphoCellsorter is a rapid way to determine if a treatment, genotype or location of cells affects morphology without going through the laborious work of 3D reconstruction. Our tool can thus be complementary to methods that provide those precise metrics.

In summary, MorphoCellSorter is an efficient and user-friendly tool for morphological analysis that operates with minimal computation time. The code is open-source on GitHub, accessible by the community (https://github.com/Pascuallab/MorphCellSorter) and improved according to the needs. It can be complemented with additional measures, such as the expression of specific markers, to investigate the relationship between morphology and gene expression. The emergence of spatial transcriptomics is a promising avenue for furthering our understanding of the link between morphology and phenotype in the microglia field in which MorphoCellSorter will be of great interest.

Availability of data and materials

The original code of MorphoCellSorter is available on GitHub (https://github.com/PascualLab/MorphocellSorter).

The datasets used during the current study are available from the corresponding author on reasonable request.

Additional information

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Funding

This study is supported by research grants from Fondation pour la recherche médicale (reference DQ20170336751) and has been developed within the BETPSY project, which is supported by a public grant overseen by the French national research agency (Agence nationale de la recherche, ANR), as part of the second “Investissements d’Avenir” program (reference ANR-18-RHUS-0012). This work was also supported by the French National Research Agency within the framework of the LABEX CORTEX ANR-11-LABX-0042.

Authors’ contributions

SB, SMD, PM, JC and OP conceptualized and designed the study. JC and OP were in charge of overall direction. The experiments were conducted by SB, MTB, LC, BP, GP CL, IH, VH, MW and SM. JH and OP secured the fundings. Coding was executed by SB, SMD, PM and JC. SB, PDG and KC participated in the analysis. All authors provided critical feedback and helped shape the research, analysis and manuscript. SB and OP drafted the manuscript, and all authors provided critical review, commentary or revision.

Acknowledgements

We acknowledge the contribution of SFR Santé Lyon-Est (UAR3453 CNRS, US7 Inserm, UCBL) facility: CIQLE (a LyMIC member), especially Bruno Chapuis and Denis Ressnikoff for their help.