Network dynamics-based subtyping of Alzheimer’s disease with microglial genetic risk factors

Abstract

Background

The potential of microglia as a target for Alzheimer’s disease (AD) treatment is promising, yet the clinical and pathological diversity within microglia, driven by genetic factors, poses a significant challenge. Subtyping AD is imperative to enable precise and effective treatment strategies. However, existing subtyping methods fail to comprehensively address the intricate complexities of AD pathogenesis, particularly concerning genetic risk factors. To address this gap, we have employed systems biology approaches for AD subtyping and identified potential therapeutic targets.

Methods

We constructed patient-specific microglial molecular regulatory network models by utilizing existing literature and single-cell RNA sequencing data. The combination of large-scale computer simulations and dynamic network analysis enabled us to subtype AD patients according to their distinct molecular regulatory mechanisms. For each identified subtype, we suggested optimal targets for effective AD treatment.

Results

To investigate heterogeneity in AD and identify potential therapeutic targets, we constructed a microglia molecular regulatory network model. The network model incorporated 20 known risk factors and crucial signaling pathways associated with microglial functionality, such as inflammation, anti-inflammation, phagocytosis, and autophagy. Probabilistic simulations with patient-specific genomic data and subsequent dynamics analysis revealed nine distinct AD subtypes characterized by core feedback mechanisms involving SPI1, CASS4, and MEF2C. Moreover, we identified PICALM, MEF2C, and LAT2 as common therapeutic targets among several subtypes. Furthermore, we clarified the reasons for the previous contradictory experimental results that suggested both the activation and inhibition of AKT or INPP5D could activate AD through dynamic analysis. This highlights the multifaceted nature of microglial network regulation.

Conclusions

These results offer a means to classify AD patients by their genetic risk factors, clarify inconsistent experimental findings, and advance the development of treatments tailored to individual genotypes for AD.

Background

Alzheimer’s disease (AD) is complex and heterogeneous, exhibiting a wide spectrum of clinical phenotypes and degrees of pathology, resulting in variability in response to treatments [1,2,3,4]. Heterogeneity in AD poses a significant challenge in developing effective therapies, as treatments that target specific molecular pathways are only effective in a subgroup of patients [5]. Precision medicine approaches, therefore, warrant attention to dealing with heterogeneity among AD patients, leading to the subtyping of AD patients [6,7,8]. However, the subtyping of AD patients has been focused on disease characteristics based on the amyloid hypothesis, which posits that the accumulation of amyloid beta plaques is the primary driver of AD pathology [9,10,11,12]. However, such an approach has limitations in explaining the inconsistencies between pathological features and clinical phenotypes observed in AD patients, which are thought to originate from the complex and nonlinear characteristics of AD pathogenesis [13, 14]. Therefore, an improved patient subtyping method that utilizes the multifactorial causes and complex interactions of AD is required to overcome this situation.

Neuroinflammation and synapse loss are recognized as key elements in AD progression among its various features. In this respect, microglia, the primary immune cells of the brain, have gained particular interest as a potential target for drug development in AD [15,16,17,18]. Because microglia play a pivotal role in AD through their inflammatory and phagocytotic functions, understanding the functions and mechanisms of microglia in AD progression is now considered a crucial task [19]. For this reason, recent studies have tried to target microglia processes, especially inflammatory processes, to attenuate AD progression [20, 21]. However, while targeting microglial pathways has shown improvement in cellular and transgenic AD model organisms with specific mutations, clinical trials have reported limited success [22, 23]. The limited efficacy of such studies stems from their focus on a single pathway or individual molecules, which may be effective only in a subset of the heterogeneous AD population [24]. It is, therefore, necessary to understand the nonlinear and complex molecular interactions within microglia in order to subtype AD patients by pathogenesis mechanisms and identify potential drug targets.

Other than the complexity of molecular interactions themselves, microglial heterogeneity can also be a large obstacle to classifying AD patients. The heterogeneity of microglia between patients is known to play a crucial role in driving diverse microglial responses in the progression of AD [25]. The largest source of microglial variability originates from genetic variability and is known to be correlated with amyloid response [26, 27]. Existing methods for assessing genetic variation employ linear approaches, neglecting the nonlinear characteristics of gene expression [13]. Consequently, while these techniques may predict the global trend in relationships between severity and each AD characteristic, they fail to uncover detailed nonlinear mechanisms underlying pathogenesis [28, 29]. Therefore, given the large influence of genetics on microglial variation and the complexity of gene expression, the integration of genetic risk factors and their nonlinear interactions is essential for subtyping AD patients [30].

In recent years, systems biology approaches have emerged as powerful tools for studying complex diseases, such as cancer and aging [31,32,33]. These approaches utilize network-based mathematical models and quantitative simulations to analyze the molecular dynamics underlying disease processes. Boolean models, which represent molecular interactions as logical operations with binary variables, are known to appropriately capture and describe biological phenomena [34, 35]. Such Boolean models have been proven to be effective in analyzing large-scale networks and explaining cellular dynamics to unveil hidden mechanisms [36,37,38]. The Boolean modeling framework can integrate genomic alterations, signaling pathways, and molecular interactions to predict how changes at the molecular level affect cellular phenotypes [39, 40]. Therefore, by employing the Boolean modeling framework, it is possible to subtype AD patients based on their genomic alterations and understand the complex molecular interactions that drive AD heterogeneity.

In this study, we constructed a comprehensive molecular regulatory network model of microglia using the Boolean modeling framework, including information on AD risk factors from scattered molecular interaction data and single-cell RNA sequencing data of microglia. The network model incorporates known major biological processes, including the inflammatory response, cytoskeleton organization, and lipid homeostasis, to explain AD-related phenotypes such as amyloid plaque formation, neuroinflammation, and phagocytosis. We integrated patient-specific mutational backgrounds into the microglia network models and simulated them until they reached stable states. Subsequently, we applied clustering algorithms to the simulated results, which revealed subtypes of AD patients based on their combination of amyloid plaque formation, inflammation, and lipid homeostasis pathways. Network analysis allowed us to identify feedback mechanisms involving SPI1, CASS4, and MEF2C, which determined AD subtypes. Furthermore, network dynamics analysis suggested PICALM, MEF2C, and LAT2 as common therapeutic targets to restore a normal phenotype in several subtypes and provided explanations for previous contradictory experimental results, which suggested both activation and inhibition of AKT or INPP5D could ameliorate AD. (Fig. 1 explains the overall methodology presented in the paper)

Fig. 1

Overall schematic of the methodology. The overall process of a microglia network-based patient classification. A microglia network including known genetic risk factors is constructed using single-cell RNA sequencing data along with literature surveys. The constructed microglia network, along with a perturbation probability obtained from SNP data, is utilized to create a subject-specific network model. Using probabilistic Boolean simulations, node activity analysis is performed for subject classification, and subtype-specific potential therapeutic targets are obtained using network simulations

In summary, we created a mathematical microglial network model to include major microglial functions. Subtyping AD patients using the constructed microglia network model allows us to uncover novel insights into the disease’s underlying heterogeneity and identify potential drug targets for personalized treatment approaches that have the potential to aid clinical decision-making.

Methods

Single-cell RNA-seq data preprocessing

We obtained public single-cell RNA expression data from the Synapse software platform from Sage Bionetworks under the project syn21438358 [41]. The dataset encompasses thirteen subjects, including those diagnosed with mild cognitive impairment or AD between Braak stages one and three.

We employed Cellranger version 5.0.2 by 10X Genomics to process the downloaded raw sequencing data, utilizing the human reference genome GRCh38 [42]. The output consisted of a cell-by-gene matrix, with each element (i, j) representing the unique molecular indices detected in cell i for gene j.

Subsequent analyses were conducted using Scanpy version 1.6.0 within the Python environment [43]. To ensure data quality, we initiated basic filtering steps. Genes detected in fewer than three cells and cells expressing fewer than 200 genes were excluded. Additionally, we filtered out potential cell doublets and apoptotic cells by removing cells with over 2,000 expressed genes or more than 10% mitochondrial reads. These filtering steps resulted in a final dataset comprising 13,324 cells, down from the initial 14,339 cells.

In order to infer core transcription factors and their interactions, the obtained filtered cells × genes matrices were concatenated into a single matrix and were used as the input for the pySCENIC package version 0.11.0 [44]. Given the stochastic nature of pySCENIC, we performed ten runs and retained transcription factor-target gene (TF-TG) pairs that were directly annotated and appeared at least eight times in these ten runs. The obtained TF-TG pairs, along with the cells by genes matrices, were utilized to quantify regulon activity using the AUCell module within the pySCENIC pipeline. Each regulon was binarized using a Gaussian mixture model to facilitate downstream analysis.

Network construction

We conducted an extensive survey of experimental and analytical data to identify genetic risk factors associated with AD in microglial cells. Our study focused on the phagocytotic/endocytotic and microglial activation pathways, which are recognized as core functions of microglia in AD [20]. We utilized core transcription factors and established connections using single-cell RNA sequencing data to establish a comprehensive network linking the identified genetic risk factors.

The initial step involved creating a gene regulatory network that incorporated the risk factors and transcription factors within one degree of separation from the identified risk factors. To ensure robust network construction, we removed self-loops to prevent the logic from being influenced by the initial state. Transcription factors that formed a weakly connected component were retained, while those transcription factors not associated with genetic risk factors and with an in-degree of zero were excluded. For each regulon, if the in-degree exceeded 12, we selected the top 12 links based on the highest average importance value obtained from the GRNBoost2 module of the pySCENIC pipeline. This selection was made due to computational constraints during the logic inference process. The remaining connections served as the foundational structure for our regulatory network.

Using the network structure and binarized regulon activities, we established Boolean logic for each transcription factor. Initially, the Boolean logic for the gene regulatory network was determined using the Quine-McCluskey algorithm, assuming that each cell was in a steady state. These inferred logics were then adjusted to create monotonic Boolean functions that optimized accuracy. In instances where the logic for a specific regulon resulted in either TRUE or FALSE values and the regulon was not a risk factor, the corresponding node was removed. Additionally, links absent in the resulting logic were removed, and isolated transcription factors were excluded from the final network.

Subsequently, we conducted extensive surveys to gather regulatory information regarding genes and proteins associated with the nodes in the gene regulatory network and the risk factors. Proteins forming complexes to activate downstream targets were connected using “AND” logic gates, while logics for risk factors obtained from different experiments or contexts was fit with “OR” logic gates. The validity of our network model was confirmed by ensuring that it accurately replicated pathological input-output relationships of microglia.

Obtaining risk alleles

To identify the risk alleles associated with AD pathogenesis within the constructed microglial network, we utilized the NHGRI-EBI-GWAS catalog [45]. A keyword search for ‘Alzheimer’ was employed to filter traits linked to the risk genes incorporated in our network. The variants within these genes were further categorized based on their statistical significance, as indicated by p-values. Our analysis concentrated on the top-ranking variants, and from the top three variants for each gene, we selected the variant with the highest risk allele frequency. This approach ensured that the chosen variants held strong potential for impacting AD risk.

To determine the number of risk alleles for each subject, we obtained whole-genome sequencing data from the Religious Orders Study and Rush Memory and Aging Project (ROSMAP), and the Mount Sinai Brain Bank (MSBB) projects. In order to obtain the effect of how risk allele affects gene expression, we examined whether these alleles led to increased or decreased gene expression. If no data was present, we identified whether the gene expression was increased or decreased in AD patients compared to their control counterparts. If a gene expression for a provided risk allele was increased in AD patients, we assumed increased gene expression would lead to increased AD risk.

Probabilistic boolean simulations

To incorporate the impact of risk alleles and their quantity within the network, we introduced a probabilistic Boolean framework capable of accommodating both gene expression upregulation and downregulation. For risk variants associated with decreased gene expression, we calculated a probability by multiplying the count of risk alleles by a predefined value (0.05). We employed this probability to switch the node state to FALSE. Conversely, for risk variants linked to increased gene expression, we computed a probability using the same multiplication approach and set the node to an ON state with the derived probability. Nodes APOE and CLU were set to a TRUE logic unless inhibited by the risk allele count, as we assumed such genes were considered external factors. In cases not falling into these categories, we updated the gene’s state according to the provided Boolean logic synchronously.

Our simulations were initiated by generating 1,000,000 random initial states, aiming to capture the diverse microglial states within a single patient. Subsequently, we simulated the network dynamics over 500 discrete time steps, with our analysis primarily focusing on the average node expression levels during the final time steps (401 to 500) to reveal stable states and patterns within the network. This probabilistic Boolean framework enabled us to comprehensively consider the intricate interplay of risk alleles and their consequent effects on gene expression, providing insights into the network’s behavior under diverse conditions.

Subject classification

Utilizing the average node expression data obtained from our probabilistic Boolean simulations, we employed Uniform Manifold Approximation and Projection (UMAP) to visualize the high-dimensional gene expression space. Subsequently, we applied Leiden clustering to partition subjects into distinct clusters based on their node expression profiles.

To unravel the node expression underpinnings of the obtained clustering, we employed the scikit-learn decision tree classifier to identify up to four risk variants that were pivotal in determining cluster assignments. The decision tree classifier was used using the assigned cluster with the Gini impurity as the splitting criterion and with a max depth of four. Such resulted in a classification accuracy of 86.9%. The three most pivotal genes in classification were selected.

The risk variants identified through the decision tree classifier were key components in constructing the expanded Boolean network. The expanded network aimed to elucidate regulatory relationships among risk genes and their connections that defined the clusters. During this process, inputs not produced by microglia cells, such as CX3CL1 or GF, were assumed to be fixed to a TRUE value. Expanded network analysis was used to explore dynamic behaviors and feedback loops involved in risk genes in AD-associated clusters. Using available clinical data, including cognitive diagnosis, Braak stage, CERAD score data, and APOE4 allele count, statistical analysis using analysis of variance was performed between clusters.

Identifying reversion targets for each subtype

In order to elucidate interventions that could restore the AD-associated phenotype to a non-AD state, we utilized a computational model for each subtype. Each subtype used 1,000,000 random initial states to assess the average expression of phenotype nodes, including phagocytosis, autophagy, pro-inflammatory cytokines (TNFα, IL1), anti-inflammatory cytokines (IL4, TGFβ), and amyloid beta tangle formation. For each cluster, we determined the most frequent allele count for each variant in both AD and non-AD subjects. The “nominal non-AD state” was defined as a state where no genetic perturbations were present, while an “AD state” for each cluster was characterized by the differences in the most frequent allele count.

To evaluate potential reversion targets, we employed deterministic Boolean simulations utilizing the Boolnet package in R [46]. We calculated an AD score to quantify the effectiveness of each control target. The AD score was defined as the sum of the absolute differences between the nominal node expression and the average node expression post-perturbations, divided by the sum of the absolute differences between the AD node expression and the nominal node expression. An AD score close to zero indicates that the control target effectively reversed the phenotype node expression to that of a non-AD state, suggesting a promising reversion candidate. In contrast, an AD score close to one signified that the control target was ineffective in restoring the phenotype to a non-AD state, reflecting a limited therapeutic potential. While an AD score higher than one may suggest an unknown third state, such a state deviates significantly from the desired nominal state and is thus irrelevant in the context of reversion.

We systematically assessed potential control targets by fixing nodes to TRUE or FALSE values and evaluated single and double targets. The targets with the lowest AD scores were identified as the most promising reversion candidates for each AD subtype and represent possible key interventions that have the potential to reverse the AD-associated phenotype towards a non-AD state.

Results

Construction of a microglia molecular regulatory network

We sought to construct a comprehensive model of molecular interactions within microglia to elucidate the intricate processes contributing to the diverse pathological traits observed in AD. To capture the essential functions of microglia implicated in AD pathogenesis, we specifically focused on the phagocytotic, inflammatory, and anti-inflammatory pathways known to be core functions of microglia. In order to systematically analyze the crosstalk between the multiple pathways, we reconstructed a Boolean model for a molecular network of microglia based on published literature and TF-TG relationships based on single-cell RNA sequencing data (See Supplementary Information 1 for comprehensive information about the nodes, links, and Boolean logical rules governing the node states.). We inferred and selected key transcriptional factors, including BHLHE41, CHD2, IKZF1, SF1, XBP1, and YY1, for inclusion in the network (Supplementary Information 2). Additionally, we incorporated known AD risk factors into the microglial network to explore their effects on microglial pathways.

The constructed microglial network comprises 63 nodes and 214 links, each node labeled by its corresponding gene name where applicable (Fig. 2). Nine input nodes (amyloid beta, growth factors, APOE, TGFβ, TNFα, CLU, CX3CL1, and interleukins 1,6 and 4,10) represent various stimuli, activating different network model pathways. Additionally, seven output nodes (TNFα, TGFβ, interleukins 1/6, 4/10, autophagy, phagocytosis, and amyloid beta plaques) represent the phenotypes and cytokines generated by the network model.

Fig. 2

Microglia network model with risk factors. The constructed microglia network model was plotted using PowerPoint software. The network contains nine input nodes and eight output nodes. The color of the node indicates the pathway of the corresponding node. The input node is indicated with a white square, while the output node is represented by a hexagon. The color or shape of the arrow indicates the direction of the interaction. The final network contains 63 nodes and 214 links. The risk factors included in the network are represented in bold

To validate the inferred network links, we first checked the correlations between the source and target genes of each link in the microglia network using the original dataset. These correlations were much more significant than those between random gene pairs (Supplementary Information 3 A, Fig. S1, S2). We then checked whether the links in the microglia network were also inferred from other independent datasets. For this, we constructed two additional networks using single-cell RNA sequencing data: one from the Human Brain Cell Atlas and another from previously published literature which presents single-nucleus RNA sequencing data from iPSC-derived microglia-like cells, using the same methods as the original microglia network [14, 47]. When we checked the overlap of the links in these networks with the microglia network, we found that almost all the links in the microglia network were present in the newly constructed networks. These results indicate that the microglia network is robust and reliable (Supplementary Information 3B, Table S1).

To validate the accuracy of our constructed network, we conducted simulations by updating random initial states using Boolean logical rules using an unperturbed network state to emulate normal aging conditions. We compared this state to a network configuration that simulates AD by manipulating amyloid beta levels. By using known experimental results, perturbation simulations on nodes were performed, and changes in the model and experimental results were compared to ensure the appropriateness and fidelity of our network model (Supplementary Information 3B, Table S2).

Effects of multiple risk factors in probabilistic boolean simulations

Standard deterministic Boolean networks use fixed Boolean functions to determine the next state of each node, in contrast to probabilistic Boolean networks, which incorporate randomness by assigning probabilities to transitions between states. Deterministic Boolean networks are limited in their use of considering risk alleles due to the fact that the count of risk alleles for each risk factor ranges from zero to two, which cannot be implemented in deterministic Boolean models. Furthermore, due to their deterministic properties in deterministic Boolean models, combinations of risk factors are abrogated. To address this challenge, we devised three illustrative toy models, each implemented with both deterministic and probabilistic Boolean models.

The first toy model illustrated the combinatorial effect of two risk genes on a downstream node (c), where the logic gate “OR” determined the node’s state based on the inputs from upstream nodes (a) and (b). In the scenario where a risk factor upregulated one of the upstream nodes, the deterministic network struggled to discriminate the presence of an additional risk factor. However, the probabilistic framework demonstrated its ability to discern another risk factor’s presence by altering the output node’s average node activity (Fig. 3A). The second toy model showed the influence of two risk genes on a downstream node (c) but employed the logic gate “AND” for node determination. Here, the deterministic network faced difficulty distinguishing the number of risk factors when a downregulating risk factor was introduced into the system. In contrast, the probabilistic Boolean network exhibited greater flexibility in accounting for the impact of risk alleles, even in the presence of downregulation (Fig. 3B). The third toy model introduced a cascading three-node network, where the signal from an upstream node (a) propagated to nodes (b) and (c). In this setup, the deterministic Boolean network overlooked the effect of the upper-risk allele, leading to both networks producing identical outcomes (Fig. 3C).

Fig. 3

Probabilistic Boolean networks are able to discriminate the types of risk factors. Deterministic Boolean simulations produce the same results, regardless of a second Boolean activating risk allele (above). Probabilistic Boolean simulations provide distinct results and distinguish the presence of the second risk factor. (A) An example of when a third node is tied with the “OR” logic. (B) An example of when the third node is tied with the “AND” logic. (C) Three cascading nodes with two risk alleles in the upstream nodes. (D) A toy network model and its deterministic/probabilistic simulation results. Each risk allele is assumed to activate or inhibit the corresponding node by 50%

In a scenario where a feedback loop is present between two nodes (a) and (b), each with a downregulating risk factor, and an upregulating risk factor, deterministic simulations fix the average node expression levels to zero or a value of one (Fig. 3D). In deterministic simulations, node (c) would be in an “OFF” state and would lead to limited signal transduction, which would be unrealistic due to risk alleles rarely having deleterious effects on the corresponding risk genes.

These findings highlighted the limitations of deterministic Boolean simulations, as they fell short in representing the complex interplay of risk alleles in AD. In contrast, the probabilistic Boolean network emerged as a more promising approach, capable of considering both the location of risk genes and the count of risk alleles, thus offering a more accurate portrayal of the intricate relationships between risk factors in AD. This distinction is crucial, especially considering that the number risk alleles is known to elevate the risk of AD, and that risk alleles rarely induce complete inhibition or activation of its corresponding risk genes. The toy model signifies the importance of adopting a probabilistic framework for a more comprehensive understanding of AD pathogenesis.

Risk variant identification and patient-specific probabilistic boolean network simulations

In our study, we obtained risk alleles for each subject and obtained perturbation probabilities. We then perturbed corresponding risk genes to create subject-specific network models, as the construction and simulation of subject-specific network models allow the identification of differences in network dynamics. These networks were then subjected to probabilistic Boolean simulations to identify the specific combinations of risk factors influencing the cytokine profiles and phenotypes of microglia. Our investigation encompassed a total of 20 risk factors, well-established for their impact on AD risk (Table 1). Information regarding these factors, including risk allele location, frequency, odds ratio, risk allele identity, and reference allele, were sourced from the NHGRI-EBI-GWAS catalog [45]. Among these factors, CLU and APOE are exogenous agents not naturally produced by microglia but were included as input nodes within the microglia network model, given their significant influence on AD risk.

Table 1 Risk genes included in the network and information used risk variants

Subsequently, we generated subject-specific networks with varying logic, incorporating data on the number of risk alleles obtained from VCF files originating from the ROSMAP and MSBB projects (Supplementary Information 3 C) [74, 75]. We then incorporated an activation or inhibition ratio of 5% multiplied by the risk allele count for risk genes for the subject-specific networks. For example, if a risk gene has one risk allele, the activation or inhibition probability would be 5%, while if a risk gene has two risk alleles, the probability would be set to 10%. The selection of a 5% probability was guided by its ability to differentiate clusters without fixing average node activity to binary values of 0 or 1. APOE, owing to its significant role in AD risk, was assigned a higher probability of 15% for its regulatory effects during probabilistic network simulations. This approach, featuring upregulation or downregulation ratios, offered greater precision, acknowledging that risk alleles often induce partial increases or decreases in the expression of risk genes.

After setting the activation or inhibition probabilities, risk genes were stochastically updated based on one of the two logics within our probabilistic Boolean modeling framework. In the case of inhibitory risk alleles, the probability of inhibition was set at 10% when two risk alleles were present. With this defined probability, the risk node had a 10% chance of being turned into a FALSE value in the subsequent update step. Otherwise, the risk node followed the original logic of the node.

Leveraging probabilistic simulations allowed us not only to account for the count of risk alleles but also to simultaneously consider multiple risk genes. We employed 1,000,000 initial states for 500 steps of synchronous updates, with the average node expression at time steps 401 to 500 selected as a suitable measure, signifying sufficient convergence of network expression levels (Supplementary Information 3D, Fig. S3). Such simulations enabled the identification of combinations of risk factors influencing microglial phenotypes and cytokine production.

Classification of subjects based on subject-specific network dynamics

The classification of subjects based on network dynamics hinges on genetic alterations, constituting the most significant source of heterogeneity within microglia, consequently impacting their transcriptional signatures. Hence, for precision medicine, it becomes imperative to perform subject classification that incorporates the genetic backgrounds of individuals, which is possible using network simulations (Fig. S4).

Our approach centered on evaluating the average activity of each node in subjects to facilitate their classification. Based on these average node activities within the network, we successfully identified nine distinct subtypes among the subjects with different clinical phenotypes (Fig. 4A and B, Table S3).

Fig. 4

Identified patient subtypes and their characteristics. (A) The result of UMAP clustering on the average node expressions, using a resolution of 0.1, and their corresponding disease status. Each point represents a single subject obtained using the average node activity post-simulation. (B) The number of APOE4 alleles and the number of AD subjects are presented for each patient subtype. (C) A heatmap of the expressions with hierarchical clustering. Each row represents the average node activity for each subject obtained via probabilistic Boolean simulations. The columns represent the nodes included in the network

These subtypes can be broadly categorized into two major groups. The first category encompasses subtypes 0, 2, 3, 5, and 8 and is characterized by a pronounced association with inflammatory signatures. Conversely, the second category comprises subtypes 1, 4, 6, and 7, displaying non-inflammatory signatures (Fig. 4C). This classification aligns with earlier findings, which also reported a similar distribution of patients within these two groupings [76]. Of particular note, the high-inflammatory group exhibited lower activity in phagocytosis nodes, higher activity in autophagy nodes, and relatively reduced levels of amyloid beta plaques. The average expression of amyloid beta plaques between the clusters correlated with Braak staging, showing that the clusters appropriately represent clinical phenotypes. (Supplementary Information 3E, Fig. S5)

Further observation of APOE4 allele count revealed that subtype 2 displayed a higher presence of APOE4 alleles, and subtype 0 exhibited a lower presence of APOE4 alleles. However, other than the two subtypes, the overall impact of the APOE4 allele on subtype classification within the microglia network was relatively minor. The cognitive diagnosis score was highest AD in subtype 2 and lowest in subtype 8. Braak staging, on the other hand, was highest in subtype 8, while others showed small differences. CERAD score was lowest in subtype 2 and high in subtypes 0 and 7. Such results show that the subtypes identified through simulations show distinct clinical phenotypes (Fig. 5). These findings are consistent with previous reports that Braak staging does not completely match with AD diagnosis [77].

Fig. 5

Clinical phenotypes for each identified patient subtype. The clinical phenotypes, including CERAD score, Braak staging, and cognitive diagnosis, are represented (* p-value < 0.05, ** p-value < 0.01, *** p-value < 0.001) for each phenotype. In the case of the APOE4 allele count, the significance levels between subtype 0 and all other subtypes are significant (p-value < 0.001) and are not represented in the diagram for simplicity

To further validate our results, we compared the subtypes identified by our method with those previously identified by Neff et al. using RNA-sequencing data [78]. Neff et al. identified five subtypes, each with distinct characteristics related to amyloid beta binding, tau-related genes, and immune pathways. We found that our network dynamics-based classification results are largely consistent with previous classification results (Supplementary Information 3 F), indicating the robustness and reliability of our subtypes based on genetic risk factors.

Network analysis using an expanded Boolean network reveals core feedbacks for each subtype

Our investigation led to the identification of core feedback mechanisms that underpin the classification of subjects within each subtype, with a particular focus on the involvement of risk genes. In this context, mutations commonly present and instrumental in classifying the subtypes were utilized to discern the feedback loops involving risk alleles using an expanded Boolean network, which integrates the regulatory rules of each node in the original network (Fig. 6A and B) [79]. Three pivotal genes -SPI1, CASS4, and MEF2C- emerged as the key determinants in subject classification.

Fig. 6

Expanded network representation of core feedbacks for distinguishing subtypes. (A) A toy Boolean model and its corresponding logic. (B) The constructed expanded network. If node 1 has an activating risk factor and node 5 has an inhibiting risk factor, the core feedback is represented in red. The dots in the network represent a composite node, in which all the upstream nodes must be in an “ON” state in order for the composite node to be turned on. Core feedbacks of (C) subtypes 0, 2, 3, (D) subtypes 5, 8, and (E) subtypes 1, 4, and 7 are presented

For subtypes 0, 2, and 3, risk alleles in the SPI1 gene were identified, while no risk variants in the CASS4 gene were observed. It is noteworthy that SPI1 is a well-known master regulator of the inflammatory response in microglia [72, 80, 81]. The presence of risk variants in the SPI1 gene played a central role in creating a positive feedback loop. This loop had cascading effects on key nodes within the microglial network, including ROS, XBP1, and NFκB. Consequently, it activated downstream inflammatory nodes, such as IL1 and TNFα, while concurrently downregulating nodes associated with anti-inflammatory processes, such as IL4 and autophagy. Importantly, the absence of risk variants in CASS4 correlated with an absence of expression increase in the CASS4 node. This lack of expression enhancement was linked to a decrease in SRC expression. Additionally, the upregulation of SPI1 further contributed to a positive feedback loop by downregulating SYK expression, leading to the inhibition of phagocytosis levels, a crucial microglial function. In subtype 3, which lacked risk variants in the ABI3 gene, a distinct positive feedback loop involving ABI3 was present, and downregulated AKT expression was present. Simultaneously, increased activity in SPI1 led to elevated CD33 and PTPN6 node expressions while reducing SYK node expression. This decline in SYK expression further contributed to the positive feedback initiated by CASS4, SRC, LAT2, PLCG, and calcium ion nodes, further decreasing phagocytosis levels (Fig. 6C).

Subtypes 5 and 8 exhibited an upregulation of CASS4 and SPI1 due to the presence of risk alleles. Additionally, PTPN6 displayed an upregulation driven by the cascade of CD33 and PILRA. SRC exhibited an upregulation influenced by both PTPN6 and CASS4. The interplay among CASS4, SRC, and PTPN6 established a positive feedback loop that led to enhanced upregulation of the phagocytosis node. Furthermore, consistent with previous observations in other subtypes, the SPI1 feedback mechanism in subtypes 5 and 8 exerted inhibitory effects on autophagy and IL4 expression. Conversely, this same feedback loop activated the expression of inflammatory markers IL1 and TNF α (Fig. 6D).

In subtypes 1, 4, 6, and 7, no risk variants in the SPI1 gene were detected. Consequently, both nodes were deactivated due to the positive feedback loop between SPI1 and NFκB. This deactivation had cascading effects on downstream nodes, leading to the deactivation of BIN1, INPP5D, IRF1, and IKZF1 nodes while activating the IRF3 node. Specifically, in subtypes 4 and 6, where MEF2C risk alleles were absent, a cooperative effect between IRF3 and MEF2C was observed. Together, they activated the TGFβ node, initiating a positive feedback loop between TGFβ and MEF2C. This led to the deactivation of downstream inflammatory nodes, such as IL1 and TNF α, while activating the anti-inflammatory node IL4. In subtypes 6 and 7, where CASS4 risk variants were present, the activation of PTPN6 triggered a positive feedback loop involving CASS4, SRC, LAT2, PLCG, ERK, CD33, and PTPN6. This loop was further strengthened by the inhibition of PI3K, caused by the activation of PTPN6. Additionally, the inactivated SPI1 node led to the deactivation of the ROS node, which, in turn, further inhibited NFκB. This series of events reinforced the positive feedback loop centered around SPI1 (Fig. 6E).

To further verify the significance of the pivotal risk genes in subtype classification, we employed other independent publicly available gene expression data obtained from 35 subjects to confirm whether the clustering using pivotal risk genes, including SPI1, CASS4, MEF2C, ABI3, and INPP5D, on RNA sequencing data shows similar phenotype patterns to our network-based subtypes [82]. We found that by using the pivotal risk genes, phenotype marker patterns for phagocytosis, autophagy, inflammation, and anti-inflammation show similar patterns to the corresponding predicted node activities. This indicates that the pivotal risk factors play an important role in subject classification (Supplementary Information 3G, Fig. S6).

Subtype classification analysis unveiled intricate feedback loops within the microglial network shaped by the combinations of specific risk variants associated with AD. These feedback loops exert influence on the inflammatory, anti-inflammatory, phagocytotic, and autophagy responses within microglia.

Identification of control targets for each subtype to restore AD to a normal phenotype

Based on the previous results, we embarked on the critical task of identifying control targets within the microglial network that could potentially be manipulated to restore a normal phenotype in the context of AD. Our approach was rooted in network perturbation analysis and the quantification of specific scores to gauge therapeutic effectiveness, which was shown in previous studies that network-based stratification is effective in predicting therapeutic responses.

To assess the impact of network perturbations, we employed a phenotype score that considered the average expression of output nodes (IL1, IL4, TGFβ, TNFα, phagocytosis, autophagy), which represent key phenotypic characteristics. The goal was to calculate a score that could indicate how close a perturbed state was to a nominal AD state (phenotype score close to 1) or to a nominal normal state (phenotype score closer to 0). An AD state was defined for each subtype using the differences in risk allele count (Supplementary Information 3 H). The phenotype score was calculated as the sum of the absolute differences between the average node expression of a nominal state, without any perturbations present, and the average node expression post-perturbations, divided by the sum of the absolute differences between the AD node expression and the nominal node expression of the output nodes as shown below, where E_normal, E_control, and E_AD represents the normal state node expression, node expression post control, and the node expression in an AD state respectively.

$$\:\text{p}\text{h}\text{e}\text{n}\text{o}\text{t}\text{y}\text{p}\text{e}\:\text{s}\text{c}\text{o}\text{r}\text{e}\:\:=\:\frac{\sum\:\left|E\_normal\:-\:E\_control\right|}{\sum\:\left|E\_normal\:-\:E\_AD\right|}$$

We selected network perturbations with the smallest phenotype score, using up to two perturbations for each subtype (Fig. 7; Table 2). The network simulations identified the inhibition of PICALM as a single target and the dual inhibition of LAT2 and PICALM for a double-target approach for the largest subtype, subtype 0. These targets were chosen based on their potential to counteract the distinct characteristics of cluster 0. Compared to control subjects, individuals with AD in cluster 0 exhibited elevated levels of phagocytosis and amyloid beta plaques, two key pathological features associated with the disease. Both single and double-target strategies reduced phagocytosis and amyloid beta plaques, effectively reverting the phenotype towards a more normal-like state (Detailed data on the other clusters can be found in Supplementary Information 3I, Table S4, S5). A literature search on the identified targets verified their potential as therapeutic targets (Table S6).

Fig. 7

AD scores and phenotype node activities for restoring a normal phenotype. (A) A bar plot showing the effectiveness of the restoration targets. A score close to 0 means that the phenotype was closer to that of a normal state. (B) The activities of the phenotype nodes for subtype 2 are represented in a radar chart. The phenotypes of an AD state (left), the effects of a single target (middle), and the effects of a double target (right) are presented, with the normal phenotype being represented in blue, and the AD phenotype and each phenotype after control being represented in red

Table 2 Optimal single and double targets to restore a normal phenotype

Furthermore, our study explored the effects of targets such as AKT and INPP5D. While these targets were not identified as the optimal reversion targets in all clusters, our findings demonstrated that their upregulation or downregulation could effectively reduce amyloid beta plaques in different genetic backgrounds, which aligns with previous research indicating that inhibition and activation of specific factors reduced amyloid beta plaques in different contexts (Supplementary Information 3 J). This underscores the complexity of microglial network regulation and the importance of considering network-based analysis when developing therapeutic strategies for AD.

Discussion

Treating AD remains a challenge despite recognizing the therapeutic potential of microglia-targeted therapies, either by restoring microglial cells or by drugs that target microglial pathways. The efficacy of these treatments differs greatly, and the underlying genetic heterogeneity among AD patients is known to play a significant role in such variability [2,3,4]. Molecular targeted therapies, however, are only known to be effective in a subset of patients or a particular genetic strain [5]. Therefore, for appropriate AD precision drugs, patients must be subtyped by the source of heterogeneity, mainly by genetic variations. Existing AD classification methods predominantly rely on phenotypic markers, such as amyloid beta plaques, tau tangles, and synapse loss, or utilize gene expression data, and they often overlook the complex genetic risk factors that underscore the disease [78, 83,84,85]. Furthermore, subtyping methods that utilize genomic data lack a sufficient explanation for the cause of the subtypes due to the use of deep neural networks, leading to difficulty in identifying therapeutic targets [86]. Therefore, it is necessary to establish a new method to subtype AD patients based on the complex mechanisms of pathogenesis using AD risk factors. To address this limitation, we have devised a probabilistic Boolean simulation method, which integrates risk factor information to categorize AD patients into distinct subtypes. Furthermore, this approach aids in the identification of potential therapeutic targets aimed at restoring normal microglial phenotypes. Our simulations have unveiled potential drug targets for subtypes, each substantiated through existing literature.

Utilizing systems biology approaches and computational simulations has proven instrumental in elucidating the intricate mechanisms underpinning complex biological phenomena and holds significant promise in patient subtyping [86]. Prior simulations for subtyping and analysis have been performed using deterministic Boolean simulations or network propagations [87, 88]. Ordinary differential equation models, however, are limited in scale of the network due to the many parameters required to be fitted and insufficient data. Conversely, Boolean models offer computational feasibility for large-scale networks. Therefore, Boolean models are effective in large-scale networks due to being computationally compliant whilst being able to sufficiently represent biological phenomena [89, 90]. Deterministic Boolean models, as shown in the result section, are incapable of simulating the number of risk factors, which are not binary values. Furthermore, deterministic Boolean models abrogate upstream risk factors due to their deterministic binary properties, leading to the inability to consider multiple risk factors simultaneously. In response, our work introduces a probabilistic Boolean network and illustrates that probabilistic Boolean networks are capable of handling the dynamics of multiple risk factors.

We evaluated subtypes using 20 risk factors related to microglia because most of the risk factors are related to microglial functions or are expressed by microglial cells [91]. We computationally identified SPI1, CASS4, and MEF2C as pivotal risk factors that distinguish the nine subtypes. SPI1, which was found to separate the nine subtypes into two large groups, is known to affect inflammatory responses, which is consistent with previously known subtyping results of AD. SPI1 is known to be a master regulator in microglial development and activation, with minor changes affecting the microglial transcriptome greatly, possibly suggesting why SPI1 is not identified in differentially expressed gene analysis [80, 81]. CASS4, while little investigated, is known to interact with FAK, affecting cellular adhesion, migration, motility, and maintaining cellular homeostasis [92]. MEF2C is a myocyte enhancer factor expressed in the nervous system and is known to prime microglia and act as an immune checkpoint, which limits the overactivation of microglia [93]. The three main identified risk factors distinguishing AD patients are highly related to microglial functions, with MEF2C and SPI1 being transcriptional factors.

Therapeutic interventions for each of the AD subtypes were also investigated. MEF2C, PICALM, and LAT2 were a few of the targets commonly identified among the several subtypes. Upregulation of MEF2C and downregulation of PICALM, identified through simulations, are also known to be possible drug targets in treating AD [94]. LAT2 is known to be highly expressed in AD and is thought to be a potential drug target [95]. Interestingly, genes such as AKT or INPP5D, while not identified as the best targets for the clusters, were identified as reversion targets when upregulated in some clusters, while a reversion target when downregulated in others. The simulations of upregulation and downregulation of targets, such as AKT, resulted in reducing amyloid beta plaques, showing that the targets depended greatly on the genetic background. These context-dependent findings are consistent with previous reports in which upregulation or downregulation of the targets in different strains resulted in distinct AD phenotype alterations [96, 97].

Our results can be leveraged to perform drug repurposing. First, our proposed subtype-specific core regulators can be examined to see if they might be targeted by existing drugs. Interestingly, the core regulators such as MEF2C, PICALM, and LAT2 are known to have the potential to modulate AD pathology and microglial functions. In addition, our microglia mathematical model can be used to quantitatively predict the impact of activating or inhibiting any of the nodes included in the model on the phenotype. By simulating the cellular response to various existing drugs with known targets, we might be able to identify the most promising drug candidates. Furthermore, by incorporating individual genetic variations for each patient into the model, we could suggest the optimal drug repurposing candidates tailored to each patient. These approaches present promising avenues for future research.

The presented method has assumptions and resulting limitations. In selecting the risk factors, only a single risk allele was selected for simplicity of computation and given identical upregulation and downregulation probabilities. Multiple risk alleles are present for a single gene, and each risk allele affects AD differently. Such factors may be considered to create a probability that more accurately reflects actual biological phenomena. The scale of the network may also be a factor that may affect the result of the classification of subjects. Other risk factors other than the 20 risk factors included in the microglia network may be included for more accurate results by expanding the network to encompass not only microglia but also other cell types. Additional analysis may be performed using matched GWAS and gene expression samples to further validate the subject classification methods.

We created a microglia network that can explain previously known experimental results and can be used to predict phenotypes by simulating utilizing the network. The advantage of using a probabilistic Boolean network is that it may be applied to other diseases where multiple risk factors play a role in disease progression. By constructing a probabilistic Boolean network and performing simulations, we demonstrated that the results could distinguish the numbers and the type of risk factors and classify subjects by their mechanisms. While our approach was applied to AD and microglia, probabilistic simulations may subtype other diseases or predict drug responses to discover other combinatorial therapeutic targets for each subject that restore the disease phenotype to a non-diseased state.

Conclusions

In this study, we introduced a method for the classification of AD patients using genetic risk factors by utilizing a microglia network. First, we constructed a network model of microglia through the integration of single-cell RNA sequencing data and a comprehensive literature search. We then employed subject-specific microglial networks derived from genetic risk allele counts and conducted probabilistic Boolean simulations. This methodological framework allowed for the identification of nine distinct AD subtypes, each characterized by core feedback mechanisms instrumental in subtype determination. We pinpointed SPI1, CASS4, and MEF2C as central risk factors contributing to the delineation of AD patient subtypes. Employing perturbation simulations tailored to each subtype, we identified potential therapeutic targets capable of reinstating AD phenotypes to a normal state. These results not only offer a method for patient subtyping based on genetic risk factors but also highlight the prospect of subtype-specific treatment strategies for AD. Further studies are required to validate the subtypes and the proposed drug targets for AD treatment.

Data Availability

The datasets analyzed in the current study include the dataset from the ROSMAP study (adknowledgeprotal.synapse.org/syn11724057) and the dataset from the MSBB study (adknowledgeprotal.synapse.org/syn11723899) which were used for subtyping AD patients. Microglia single-cell RNA sequencing datasets (adknowledgeprotal.synapse.org /syn21438358) were used for creating the microglia network model. The code used for probabilistic Boolean simulations is available on GitHub at https://github.com/jhchoii/Microglia_PBN.