Answer ALS is the largest collaborative effort in ALS bringing together multiple research organizations and key opinion leaders. Over 800 ALS patients and 100 healthy controls from 8 neuromuscular clinics distributed across the United States were enrolled in this project. A blood sample was collected at the first visit of each participant and iPSC lines were generated from peripheral blood mononuclear cells extracted from whole blood via an episomal iPSC reprogramming system. The consortium generated multi-omics data comprising of genomic, epigenomic, transcriptomic, proteomic, laboratory test, medical records, and other data (Baxi et al., 2022). We used transcriptomic records within the files named bulbar_vs_limb.csv and ctrl_vs_cas.csv which are currently being expanded for future research and competitions. These files were available on Kaggle for academia and industry. The former data file is meant to differentiate between how ALS initializes, specifically in the bulbar region or limbs, allowing our system to extract key sets of genes that are active in different patient subpopulations. The latter data file was used to differentiate biological mechanisms that play a role in ALS in general, and to generate genetic hypotheses about ALS subpopulations. The data used in the preparation of this article were obtained from the Answer ALS Data Portal (AALS-01184). For up-to-date information on the study and access to the data please visit https://www.answerals.org/.

An ML playground environment called NetraAI (Qorri et al., 2020; Choi et al., 2021; Cook et al., 2023) was made available to scientists at the Gladstone Institute. This allowed medical experts to interact with the ML-generated hypotheses to evaluate the findings and examine the etiological factors that were being suggested. Here, we bridged the gap that exists between advanced ML techniques and human medical expertise through augmented intelligence (Crigger et al., 2022). The methods used for the generation of the hypotheses that led to the target classes described in this paper consisted of ML methods paired with systems biology methods. In this context, we refer to ML-generated hypotheses as proposed insights about a patient subpopulation that satisfy the following criteria:

The insight must be about a specific subset of samples that the AI finds and include a multi-factor signature that pertains to this subpopulation.

An important issue is the small number of samples within the dataset used, as we did not augment our process with other data such as literature or other genetic datasets. Our process is based on authentic limitations that exist in rare disease clinical trials, which begins with inherently small sample sizes of patients. For this reason, we built an ML pipeline using methods suitable for smaller sample sizes. By allowing the algorithms to segment the patient samples into clusters of varying confidence, and extracting precisely what factors are driving each cluster, we have a set of hypotheses that can be tested statistically and by human ALS experts. Smaller datasets do not have the sample size to accurately represent the variety of manifestations of ALS, but the sample we had access to did provide insights into statistically significant patient subpopulations. The novelty of our approach stems from the following insights:

These subpopulations were then used to extract features that were supported through significance testing and expert validation. These were then used to seed biological network analyses and hypothesis generation. This is an example of augmented intelligence, where ML methods are used to enhance human expertise, especially when datasets are limited in sample size. This process was implemented as follows:

The ML methodology utilized is outlined in Figure 1 and has previously described in more detail (Qorri et al., 2020; Choi et al., 2021; Cook et al., 2023). This was the methodology used to segment the patient population before applying the biological network methods described above.

For a foundational understanding of the data’s structures and to facilitate feature reduction, a series of methods and parameters were adopted. During data preprocessing, features were centered by subtracting their respective means. Recognizing the varied feature scales, the data underwent standardization to ensure every feature converged to a mean of 0 with a standard variance of 1. When implementing the PCA, we opted for the “full” solver, a choice influenced by the manageable size of our dataset which promised a thorough decomposition. To zero in on the optimal components, a significant focus was placed on the cumulative explained variance, ensuring our emphasis was on principal components accounting for 95% of the total variance. This approach was further cross-referenced by inspecting the “elbow” of the scree plot. The significance of features was gauged through their loading values, where features with pronounced absolute values were considered for the selection process. Furthermore, these features from PCA loadings were assessed against our domain expertise. This ensured that the pruned feature set was not just technically sound but also contextually relevant, particularly in the lens of potential ALS drug targets. Before embarking on these steps, multicollinearity among features was scrutinized using the variance inflation factor (VIF). Features breaching a VIF of 10 were given a closer look. With the dataset’s size being on the smaller side, outliers posed a risk of disproportionate influence. To counteract this, data distributions were visually examined and complemented with statistical methods geared toward outlier identification and assessment.

Our study also made use of the Random Forest method, an ensemble learning technique used for its capabilities in both classification and regression tasks. By leveraging a collection of decision trees, each being trained on a randomized assortment of data subsets and features, an aggregate predictive outcome was pursued. The primary intent here was the validation of features unearthed using our unique techniques. The dataset was strategically bifurcated, earmarking 80% for training purposes and the remaining for testing. Stratified sampling was integral in this division, a necessity arising from the class imbalance observed in our target variable. Utilizing the Scikit-learn library available in Python, we initialized the Random Forest with parameters such as 500 trees, the criterion set as “gini,” max depth restricted to 30, min_samples_split and min_samples_leaf defined as 5 and 2, respectively, and finally, a consistent random_state of 42. Post training, the Gini importance was extracted, which subsequently played a pivotal role in ranking features. A predetermined threshold was set at 0.005 for feature importance, selecting only those that surpassed this benchmark. Their inclusion was further bolstered by an out-of-bag (OOB) error measuring 0.03. For evaluations, a fresh Random Forest model was trained using the cherry-picked features, which was then validated against the testing subset. Finally, grid search was utilized specifically for hyperparameter fine-tuning, resulting in optimal parameters of n_estimators at 550 and max_depth solidified at 32.

In our approach with t-SNE, we settled on settings such as a perplexity of 30. This was largely due to its alignment with smaller datasets, effortlessly balancing between local and global structures. Accompanying parameters included a learning rate of 200, capped iterations at 1000, early exaggeration of 12, a balancing angle of 0.5, and a swift PCA-based initialization for the sake of faster convergence. Additionally, the metric was strictly defined as “euclidean.” We chose “exact” for the method parameter, offering an advantage over the Barnes-Hut approximation for petite datasets, all while minimizing complexity.

HDBSCAN clustering was configured as follows: The Minimum Cluster Size was fixed at 5, with the Minimum Samples mirroring this value by default. The Cluster Selection Method was distinctly marked as “eom” or Excess of Mass. In this phase, the Allow Single Cluster option was purposefully deactivated. Alpha was precisely set at 1.0, keeping avenues open to experiment with elevated values. The metric was once again aligned with the previous selection of “euclidean,” and Core Distance was singularly set at 1 to bolster computation times.

Lastly UMAP was used to decipher the intricate interrelations among patients, echoing discoveries from our in-house methods. Crucial parameters here were the n_neighbors fixed at 15, min_dist tailored to 0.1, the metric used was “euclidean,” a spread adjusted to 1.0, and the min_dist_fraction set at 0.1 for this study.

This comprehensive approach, underlined by these carefully chosen parameters, was our roadmap to robust, interpretable results, all the while side-stepping pitfalls like overfitting and computational lags.

TargetMine, an Open Source and peer reviewed tool that uses known genetic relationships with disease, biological pathway data, and current drug information was used to provide confidence levels for the targets we discovered with our NetraAI system (Chen et al., 2022). We compiled a comprehensive list of genes, all of which are included in this study. This list was formatted into a comma-separated value (CSV) file for computational analysis. The dataset was uploaded to the TargetMine platform, where we specifically selected Homo sapiens as the reference organism. In the “Analyse Data” tab, we initiated the analysis procedure, where it was imperative to rectify the nomenclature of several genes to ensure system recognition. Following the successful recognition of all the genes, we proceeded with the detailed analysis. TargetMine generated a downloadable report, of which the disease pathway enrichment section was of particular interest to our study. This section provided the statistical significance measures that underpinned our findings and facilitated the stratification of our target genes based on confidence levels and putative functionalities. All these data including pathway enrichment provide the data to derive significance values for the targets discovered by our process.

Several studies have attempted to identify key players in ALS pathology with hopes of elucidating relevant drug targets for this fatal disease (Batra et al., 2019; Hedl et al., 2019; Nowicka et al., 2019; Wu et al., 2019). However, many identified targets relate to mitochondrial dysfunction, protein aggregation, RNA processing, axonal transport, oxidative stress, apoptosis, SOD1, phosphorylation, and the neuromuscular junction (Batra et al., 2019). Most methods to extract these targets are based on symptoms and the mechanisms of disease development and progression; however, due to the heterogeneity of the disease, it is important to identify key players that can be druggable in specific subsets of ALS patients. Several ML approaches have identified key genetic targets, and using NetraAI, we were able to verify several of the same gene targets that have been recently reported (Table 1) as well, we identified several genes that belong to the same gene family as those previously reported (Table 2; Pun et al., 2022). The functions reported in Tables 1, 2 are based on the protein family function as well as supporting literature that discusses a proposed mechanism or function. Within Table 1, DNM3TA, ERN1, HSPD1, PPIA, VCP, MAP3K5, MAKPK1, NOS1, PTK2, PTPRC, and RARA were previously identified high-confidence therapeutic targets from iPSC-differentiated motor neurons (diMN)-derived and CNS data that belonged to the druggable classes defined by PandaOmics, with supportive evidence on their ALS or neurodegeneration, and ranked as the top-50 targets in at least one of the meta-analyses (Pun et al., 2022). In contrast, PPP3CB, was identified as a novel therapeutic target in the previous reported findings (Pun et al., 2022). The findings presented in Table 2 represent gene targets that belong to the same protein family as other targets identified by the PandaOmics study.

In addition to the drug targets shown in Table 1, which have already been previously reported and validated, as well as the targets shown in Table 2, which belong to the same gene family as those previously reported (Pun et al., 2022), NetraAI was able to uncover several targets that may shed light into ALS pathophysiology and treatment efforts. Interestingly, these targets can be grouped into a collection or family, called “target classes,” that align to a unique characteristic related to ALS (Figure 2). The target classes discussed here include inflammation, epigenetic, heat shock, neuromuscular junction, autophagy, apoptosis, axonal transport, and excitotoxicity.

These targets are not exhaustive, they represent select target classes that have the potential to play a role in ALS that warrant further investigation. Collectively, these target classes suggest that the simultaneous targeting of several key hallmarks of ALS with combination targeted therapy may have the potential to slow progression, with an enhanced possibility of maintaining and sustaining an improved quality of life (QoL) for certain ALS patients.

Neuroinflammation is suggested to begin in early ALS pathogenesis, with nervous and peripheral immune systems being impacted (McCauley and Baloh, 2018). Interestingly, we were able to distinguish an inflammation target class involving TNFα (Table 3). Given the role of TNFα in immune and inflammatory activity, this is not surprising, considering an innate immune response is characteristic of neurodegenerative diseases like ALS (McCauley and Baloh, 2018). However, the role of TNFα and its receptors TNFR1 and TNFR2 are controversial, with both protective and detrimental effects being reported (Guidotti et al., 2021). Considering that neuroinflammation is a complex and atypical inflammatory process that is meant to protect the central nervous system from injury, in ALS, chronic neuroinflammation can lead to dysregulation that contributes to neurodegeneration. It is now thought that neuroinflammation has dual function, contributing to neuroprotection and possibly leading to neurotoxicity (Tortarolo et al., 2017; Guidotti et al., 2021).

Epigenetic hallmarks have been linked with ALS, specifically with histone deacetylases (HDACs) and their inhibitors, highlighting a potential therapeutic avenue for ALS patients (Klingl et al., 2021). Using the same patient dataset, we also discovered a set of candidate genes that indicate potential HDAC dysregulation and methylation (Table 4). The genes in this target class encode numerous proteins associated with DNA binding and transcription factors, particularly histones and nucleosomes. Although HDAC is a known target for several disease states, including ALS, several HDAC inhibitors currently available have a host of toxic side effects and warrant further investigation to target specific HDACs in specific patient subgroups (Janssen et al., 2010). Collectively, these results highlight the role of epigenetic regulation in ALS pathophysiology.

In ALS, motor neurons have a deficit in the ability to activate the heat shock response (HSR) and do not upregulate the expression of heat shock proteins (Hsps) which are inhibitors of apoptosis and exert an anti-inflammatory response in glia (Apolloni et al., 2019). Here, we were able to uncover a heat shock target class, where the proteins encoded by the genes of interest are primarily associated with protein transport, such as, dynein, actin, and microtubules (Table 5). Evidently, ALS is driven by a collection of genes, with cases being highly heterogeneous; however, protein aggregates in the brain and spinal cord that are positive for SOD1, TDP-43, or OPTN are present in nearly all ALS patients. Under normal physiological conditions, these protein aggregates are prevented and cleared by Hsps, providing further evidence that ALS motor neurons have an impaired ability to induce the HSR (Seminary et al., 2018).

In the context of ALS, distal axonopathy is a central hypothesis in the early stages of the disease where pathological changes occur at the neuromuscular junction (NMJ). Acetylcholinesterase (AChE) plays a crucial role in nerve-muscle contact, facilitation of neurite outgrowth, and NMJ formation and survival. Interestingly, ALS patients are characterized by abnormal AchE content in plasma, which may reflect neuromuscular disruption (Campanari et al., 2016). Here, we found HSPG2 to characterize the neuromuscular junction target class (Table 6). Interestingly, a research paper reported on HSPG2, among others, as a novel candidate mediator for disease progression. HSPG2 plays a role in immunological and inflammatory disease, neurological disease, and skeletal and muscular disorders (Morello et al., 2019).

Similar to the epigenetic target class, the accumulation of protein aggregates is proposed to disrupt cellular processes that ultimately result in neurodegeneration. Evidently, this protein aggregation in neurons is a hallmark of ALS and may be due to defects in autophagy (Ramesh and Pandey, 2017; Amin et al., 2020). Here, in the autophagy target class, we uncovered several genes implicated in the cellular processes regulating autophagy (Table 7). Autophagy is responsible for maintaining cellular and protein homeostasis in response to nutrient depletion or organelle damage (Ramesh and Pandey, 2017). However, it is still unknown whether activation or inhibition of autophagy would be most effective in the treatment of ALS (Nguyen et al., 2019). Interestingly, SOD1 is a frequent ALS mutation and it is expected that aggregation of mutant SOD1 (mSOD1) is a crucial event in ALS pathogenesis, and dysregulation of autophagy has been linked to SOD1 aggregates in motor neurons (Nguyen et al., 2019). This highlights the need to study further and identify therapeutic agents that target the clearance of these protein aggregates.

In ALS, there is evidence of apoptosis through DNA fragmentation, caspase-9 activation, BAX overexpression, and reduced Bcl-2 expression (Erekat, 2022). Interestingly, mSOD1 induces apoptosis via cytochrome c release and Bcl-2 degradation (Erekat, 2022). As a result, treatments targeting apoptosis can be helpful in rescuing neurons from cell death. In the apoptosis target class, several caspases as well as apoptotic mediators were identified (Table 8). Of note is that of the caspases identified, based on their mechanism of action and their position in the apoptotic signaling pathways, apoptotic caspases can be initiatory caspases (caspase 2, 8, 9, and 10) and executioner or effector caspases (caspases 3, 6, and 7) (Erekat, 2022). As a result, similar to autophagy, whether promoting or inhibiting critical caspases involved in apoptosis presents as a therapeutic approach for ALS patients.

Neurons have long axonal projections that rely on cytoskeletal integrity to maintain axonal stability, transport, and signaling (Theunissen et al., 2021). In ALS there is selective, early degeneration of motor neurons in the brain and spinal cord. Related to this, we identified a target class characterized by several genes that play a role in microtubule cytoskeletal organization (Table 9). Disrupted transport mechanisms can affect mitochondrial metabolism and degeneration, protein degradation, and RNA transport, collectively resulting in motor neuron death (Le Gall et al., 2020). Furthermore, within this target class, we identified TARDBP and RPA1 which have been implicated in ER-Golgi transport dysfunction that is associated with ALS (Soo et al., 2015). It is important to note that in this target class, we identify two HDACs, and conversely, in the heat shock target class, we identified dynactin. This observation demonstrates that ALS pathophysiology is characterized by overlapping systems and is heterogeneous (Le Gall et al., 2020). It should be emphasized that even though these gene candidates are organized under specific categories, the manifestation of the disorder, and the potential treatments, all depend on the fact that the corresponding proteins, and higher-order systems, interact with each other. Hence, these findings should not be considered isolated processes but parts of an emergent system.

Finally, we extracted a collection of genes implicated in excitotoxicity (Table 10). Excitotoxicity is a phenomenon that describes the toxic actions of excitatory neurotransmitters where prolonged activation starts a cascade of neurotoxicity that ultimately leads to the loss of neuronal function and cell death (Armada-Moreira et al., 2020). Importantly, excitotoxicity can both contribute to as well as be a result of other deregulations, including mitochondrial dysfunction, neuronal damage, and oxidative stress (de Marco et al., 2022). Similar to other target classes, there is evidence that dysregulation of mitochondrial calcium handling plays a role in excitotoxicity (Verma et al., 2022).

Utilizing an Open Source bioinformatics tool, TargetMine, we evaluated the confidence in the drug targets in the manuscript thus far (Chen et al., 2022). Adding our target genes to TargetMine we were provided with 11 overarching pathway categories, each with varying levels of confidence (Table 11). Of the 86 targets, 36 were associated with pathways of neurodegeneration including ALS, with a high level of confidence (3.4×10^-16). Interestingly, 30 targets were also associated with SARS-CoV infection and interferon signaling, 35 targets were associated with RHO GTPase effectors, nuclear receptor signaling, chromatin modifying enzymes and viral carcinogenesis, 61 targets were associated with nervous system development, and 43 targets were associated with homeostasis and the neuronal system, all with high levels of confidence. All of the targets identified to be associated with cell cycle, transcriptional dysregulation in cancer, organelle biogenesis and maintenance, carboxyterminal post-translational modification of tubulin, bacterial infection pathways, and autophagy, which, despite having a lower confidence level, highlight that ALS may be a complex disorder. However, an alternative explanation is that there is a historical bias toward favored pathways and that genes are inherently promiscuous, making our molecular machinery highly connected. The output of the TargetMine software is included as a Supplementary file, which includes one table outlining the statistical significance of the pathways enriched for and the other with the pathways and genes themselves.

Utilizing a dataset consisting of 31 bulbar onset and 85 limb onset ALS patients, we identified distinct subpopulations, each defined by a specific set of driving genes (Figure 3). A subpopulation of 13 limb onset ALS patients was identified to be characterized by an elevated expression of IL200RA and LRRC23 (Loop 1). Even further, we identified a distinct subpopulation of 11 bulbar onset ALS patients (Loop 2) that was characterized by a decreased expression of TBC1D20, ALG3P1, CROCC2, AC109439.1, FAM151A, and NKX2101-AS1, and an elevated expression of TMEM14A. The remaining limb onset patients, which comprised the majority of the dataset, were characterized by expression patterns opposite to the bulbar subpopulation – specifically increased expression of TBC1D20, ALG3P1, CROCC2, AC109439.1, FAM151A, and NKX2101-AS1, and decreased expression of TMEM14A. These findings indicate that specific genetic factors may accurately delineate novel subtypes of bulbar and limb-initiated ALS. Unraveling these subpopulations has significant implications for clinical trials, as it can unveil alternative etiological subtypes that might respond more favorably to particular therapeutic interventions. A gene interaction network constructed of TMEM14A and FAM151A, revealed nearest neighbor connections to RAB1, RAB2, and TDP-43 (TARDBP in the gene interaction figure), suggesting the identification of a more aggressive ALS subpopulation within the bulbar onset patients (Figure 4).

We adopted z-score normalization prior to generating the heatmap (Figure 5) facilitated by the Seaborn Python library (Waskom, 2021). It is evident that certain genes distinctly differentiate the samples across respective classes. However, a limitation of this visual representation is the inability to distinctly highlight the subpopulations present within the heterogeneous sample group. This distinction emerges prominently through ML applications, where synergistic effects arise from integrating multiple variables concurrently. Nonetheless, the distinctiveness of several genes can be ascertained by contrasting the intensities above and below the demarcating black bar in Figure 3. The heatmap corroborates the highlighted trends explained in Figure 3, specifically that TMEM14A is upregulated for bulbar-initiated samples, while TBC1D20, ALG3P1, CROCC2, AC109439.1, FAM151A, NKX2.1.AS1 are all downregulated. IL20RA and LRCC23 are upregulated for limb-initiated samples, especially for the 13 samples represented in Loop 1 of Figure 3.

We employed two classifiers, namely Random Forest and Gradient Boosted Trees, to assess their performance using a leave-out cross-validation approach. The Gradient Boosted Trees exhibited an accuracy of 70.7% in 10-fold cross-validation and 73.9% in 5-fold cross-validation, while the Random Forest classifier performed slightly better with accuracies of 74.1 and 75% in the respective cross-validation schemes. These results suggest the presence of discernible patterns within the data.

To validate the robustness of these subtype discoveries alongside the previously mentioned driving transcriptomic factors, we constructed a new dataset comprising only these relevant variables and re-evaluated the classifiers using leave-out cross-validation. Notably, the use of this reduced dataset led to enhanced model accuracy. For instance, complex models like Random Forest yielded accuracies exceeding 80% in both the 10-fold and 5-fold cross-validation iterations. Most notably, simpler models like logistic regression, which initially exhibited poor performance with an accuracy of approximately 65%, now generated stable models with an impressive accuracy of approximately 84% for both 10-fold and 5-fold cross-validations.

These findings highlight the utility of our approach in identifying subpopulations and driving transcriptomic factors, which can be further scrutinized through bioinformatics analyses. The improved accuracy of the models underscores the importance of considering these factors when characterizing ALS subtypes and devising tailored therapeutic strategies.

These targets were discovered after allowing ML to generate hypotheses about important genetic variables using the knowledge of protein–protein interactions and co-expression to extend our search. Protein interaction networks represent a rich source of data for understanding complex biological systems and deriving potential drug targets. These networks represent nodes and their interactions as edges, forming a complex graph that can be analyzed using various network analysis techniques.