The mechanistic functional landscape of retinitis pigmentosa: a machine learning-driven approach to therapeutic target discovery

Abstract

Background

Retinitis pigmentosa is the prevailing genetic cause of blindness in developed nations with no effective treatments. In the pursuit of unraveling the intricate dynamics underlying this complex disease, mechanistic models emerge as a tool of proven efficiency rooted in systems biology, to elucidate the interplay between RP genes and their mechanisms. The integration of mechanistic models and drug-target interactions under the umbrella of machine learning methodologies provides a multifaceted approach that can boost the discovery of novel therapeutic targets, facilitating further drug repurposing in RP.

Methods

By mapping Retinitis Pigmentosa-related genes (obtained from Orphanet, OMIM and HPO databases) onto KEGG signaling pathways, a collection of signaling functional circuits encompassing Retinitis Pigmentosa molecular mechanisms was defined. Next, a mechanistic model of the so-defined disease map, where the effects of interventions can be simulated, was built. Then, an explainable multi-output random forest regressor was trained using normal tissue transcriptomic data to learn causal connections between targets of approved drugs from DrugBank and the functional circuits of the mechanistic disease map. Selected target genes involvement were validated on rd10 mice, a murine model of Retinitis Pigmentosa.

Results

A mechanistic functional map of Retinitis Pigmentosa was constructed resulting in 226 functional circuits belonging to 40 KEGG signaling pathways. The method predicted 109 targets of approved drugs in use with a potential effect over circuits corresponding to nine hallmarks identified. Five of those targets were selected and experimentally validated in rd10 mice: Gabre, Gabra1 (GABARα1 protein), Slc12a5 (KCC2 protein), Grin1 (NR1 protein) and Glr2a. As a result, we provide a resource to evaluate the potential impact of drug target genes in Retinitis Pigmentosa.

Conclusions

The possibility of building actionable disease models in combination with machine learning algorithms to learn causal drug-disease interactions opens new avenues for boosting drug discovery. Such mechanistically-based hypotheses can guide and accelerate the experimental validations prioritizing drug target candidates. In this work, a mechanistic model describing the functional disease map of Retinitis Pigmentosa was developed, identifying five promising therapeutic candidates targeted by approved drug. Further experimental validation will demonstrate the efficiency of this approach for a systematic application to other rare diseases.

Background

Retinitis pigmentosa (RP) is a rare hereditary disease that mainly affects the photoreceptors in the retina, resulting in progressive vision loss. RP is typically diagnosed in childhood and rather than a single entity, RP is better described by a group of disorders characterized by night blindness, narrowing of the visual field, and ultimately, legal blindness [1]. Due to RP’s relatively high prevalence among rare diseases, affecting more than 1.5 million patients worldwide (from 1/2500 to 1/7000) [2], it is considered the most common inherited disease of the retina [3]. It is a complex disorder with multiple genetic causes and diverse clinical manifestations, counting mutations in over a hundred genes known to lead to the development of the disease [4]. However, the genetic cause of RP remains unknown in 50% of the patients [5]. In addition, current therapeutic options for RP are limited, with only a few treatments available that can slow down the progression of the disease or improve the quality of life of patients.

Increasing the knowledge of the molecular basis of the disease is a crucial step in the search for new therapeutic options to treat RP. The exponential growth of genomic data, fueled by sequencing technologies [6], has revolutionized the field of rare diseases with foreseeable achievements such as being able to diagnose all Mendelian diseases in the near future [7]. The knowledge generated in the last decade has expanded by more than two orders of magnitude the number of variants with known phenotypic effects [8]. However, the study of diseases has undergone a conceptual shift, moving from a focus on individual genes to exploring how they interact within cells [9]. This shift has led to the use of network-based models that rely on accumulated biological knowledge to identify causal relationships between genes, as defined by biological pathways [10]. The so-called “mechanistic models of human diseases” were developed by combining these novel approaches with biological human data.

Mechanistic models can be defined as mathematical representations of biological systems that aim to capture the underlying mechanisms that drive the behavior of these systems. Interestingly, these models can provide a causal link between gene expression and functional cell behavior [11]. However, the real potential of mechanistic models relies on their capability to reduce the complexity of the data while reinforcing interpretative power, unveiling specific interactions and functional outcomes [10]. Indeed, mechanistic models have helped to understand the disease mechanisms behind different cancers [12,13,14], rare diseases [15], complex diseases such as diabetes [15], or obesity [16], as well as the mechanisms of action of drugs [17], including drug repurposing [18] or gender-specific effects of drugs in cancer [19]. A relevant property of mechanistic models is that they model causality in a quite accurate manner [9]. Therefore, they can predict the potential consequences of gene perturbations (e.g. knock-out, drug inhibition effect, etc.) over pathway activities, as well as the downstream functional consequences in the cell [10, 20]. They have been used to predict candidate essential genes in cancer cell lines, whose further inhibition validated the predictions [12]. Repurposable drug candidates were proposed for Fanconi Anemia [17] which were further experimentally validated [21], and for COVID-19 [18], also validated using a cohort of around 16,000 patients of the Andalusian healthcare database [22, 23]. These computational methods represent a leap in the generation and testing of hypotheses guided by a specific functional rationale, especially useful in scenarios with limited knowledge like rare diseases such as RP.

The development of novel therapies is a slow and costly process, particularly in rare diseases, where orphan drug research accounts for limited resources. Drug repurposing, the process of discovering new therapeutic uses for existing drugs [24], emerges as a potential solution to this problem in this context [25], reducing the time and expenses involved in drug development [26].

In this work, we propose a data-driven approach to construct actionable models that describe the disease mechanism of RP. Moreover, by making use of explainable machine learning models, trained with publicly available gene expression data, we have estimated the impact that targets of drugs already in use have in the signaling mechanisms of RP, allowing us to infer and select the most suitable therapeutic target genes that will lead us to drug candidates for its repurposing in RP. Indeed, the role in RP progression of six of the proposed targets has been experimentally validated. The comprehensive analysis of the results enabled us to identify drug targets that hold significant potential in impacting disease mechanisms able to revert cell functional processes to a healthy state.

Methods

The overall methodology developed for RP modeling is divided into three main stages depicted in Fig. 1. Briefly, the process entails collecting data on RP-related mechanisms from existing open-access databases (DBs). We ensured the precision and fidelity of RP-phenotypes and RP-gene associations by using standardized definitions from OrphaNet [27], OMIM [28], and the Human Phenotype Ontology DBs [29]. The information retrieved, together with the accumulated knowledge on functional gene interactions provided by the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways DB [30], is used to produce a map of functional interactions among RP-associated genes and genes upstream and downstream in the pathways, thus defining a set of sub-pathways, the so-called signaling circuits, which can be considered a reasonable representation of the RP disease map [17, 18]. This RP mechanistic map can be used to build a signal transduction model that allows the estimation of RP-related cell functions from gene activities measured as gene expression levels [12]. The model of the RP mechanistic map is then integrated into an explainable machine learning (ML) model, which infers potentially causal relationships between genes of interest, in our case known drug targets (KDTs) extracted from the DrugBank DB, and the RP-related functional activities as described in the RP mechanistic map. The ML methodology prioritizes KDTs with a high predictive score, degranulating the results into the specific influences that each KDT has over the different parts of the mechanistic map. The results highlight KDTs that are predicted as relevant in regulating the functional landscape of RP mechanistic map and could potentially be further studied as promising therapeutic targets.

Fig. 1

Diagram of the steps involved in the machine-learning-based drug repositioning modelization of Retinitis Pigmentosa (RP). RP Disease Map construction (yellow block); machine learning modelization of drug-targets over the RP map (blue block); target prioritization for experimental validation (pink block)

Gene expression data

Transcriptomic data from the GTEx (Genotype-Tissue Expression) Portal was used to obtain meaningful interactions among proteins and the role of these proteins in the modulation of the signaling pathways. Data from 54 non-diseased tissue sites across 948 individuals, with a total of 17,382 samples over 22,000 gene expression values were downloaded from the GTEx Portal (GTEx Analysis Release V8; dbGaP Accession phs000424.v8.p2) [31]. Gene expression data from the RNA-seq GTEx dataset was normalized by applying the Trimmed mean of M values (TMM) normalization using the edgeR R/Bioconductor package [32].

Data on drug targets

To investigate the potential effect druggable molecules have over our RP-mechanistic model, the DrugBank DB release 5.1.8 [33] was used. Drug-target interactions involving human protein-target of approved drugs with a known pharmacological action were extracted. Only KDTs with gene expression measurements in the GTEx V8 dataset [31] were taken into account.

Defining the core knowledge for the mechanistic map of retinitis pigmentosa

We define the mechanistic disease map as the collection of functional and directed interaction among genes, gathered in the physiological signaling pathways, that best resembles the disease mechanism. Therefore, all the knowledge available on RP-related genes must be collected to be set in the context of functional interactions provided by the KEGG pathways. To construct the RP mechanistic map, as a first step, we extracted the genes with mutations known to be related to the development of RP from OrphaNet [27] (ORPHA:791) and OMIM [28] databases.

Expanding the mechanistic map of retinitis pigmentosa

To expand the core map, other genes linked to RP phenotypes that could be playing a key role in RP progression mechanisms were added from the Human Phenotype Ontology (HPO) database [29]. To do so, we extracted the phenotype-disease and phenotype-gene interactions from HPO phenotypic abnormality sub-ontology (hp.obo release 2019-09-06). A total of 8822 HPO-terms tagging a total of 5440 genes were assessed in this study. Specificity levels were assigned to the HPO-terms based on the hierarchical properties of the DB. Since the HPO ontologies have a hierarchical structure, the specificity level was calculated as the number of ancestor nodes an HPO-term has. Low-specificity HPO-terms can be very broad with many related processes and genes while high-specificity HPO-terms can describe very isolated ones with few connections. Only HPO-terms with a specificity level between 7 and 15, with traceable author statement (TAS) evidence sources were taken into account.

Defining the mechanistic map of retinitis pigmentosa

This set of RP-related genes (extracted from both HPO and OrphaNet/OMIM databases) was mapped to the KEGG signaling pathways using the Hipathia R Bioconductor package [34]. Then, the receptor-effector sub-pathways, also known as circuits, containing the RP-genes were extracted. The collection of circuits containing the selected RP-related genes is what we define as the RP Mechanistic Map, an actionable in-silico model where simulations of different scenarios can be performed. Those pathways with a high number of disconnected nodes or which circuits have a maximum number of nodes of three were ruled out to minimise the potential errors in signal propagation algorithm estimation of Hipathia. The complete RP mechanistic map model can be explored at: http://hipathia.babelomics.org/RP_Mechanistic_Map/ [35].

Modeling of the mechanistic map of retinitis pigmentosa

Mechanistic models provide a powerful framework for understanding the complex relationships between different components within a biological system. In this context, we used the mechanistic pathway analysis (MPA) tool Hipathia to model RP functional landscape. Hipathia is based on a recursive signal propagation algorithm that considers biological pathways as a collection of “signaling circuits”, graphs representing the interaction between proteins where the signal propagates from receptor-proteins to effector-proteins in a single path. Hipathia uses gene expression values as proxies of the level of activation of the corresponding protein node in each circuit. Taking into account the inferred protein activity and the interactions between the proteins (activations or inhibitions) defined in the pathway, the level of activity of each circuit is estimated from gene expression values using the following recursive equation.

$$S_{n} = v_{n} \cdot \left( {1 - \prod\limits_{{S_{a} \in A}} {\left( {1 - S_{a} } \right)} } \right) \cdot \prod\limits_{{S_{i} \in I}} {\left( {1 - S_{i} } \right)}$$

S_n being the signal intensity for the current node n; v_n being its normalized gene expression value; A being the set of activating signals (s_a) arriving at the current node from activation edges; I being the set of inhibiting signals (s_i) coming to the node from inhibition edges. A detailed explanation of the method can be found in Hidalgo and cols 2017 [12].

The Hipathia method, as implemented in the Hipathia R/Bioconductor package [36], was used to estimate signaling circuit activities within the 79 physiological signaling KEGG pathways from the corresponding GTEx normalized gene expression values (rescaled to [0, 1] as recommended by the package documentation).

Explainable machine learning model

This method employs a Multi-output Random Forest [37] and SHapley Additive exPlanations (SHAP) [38], tailored for tree-based algorithms [39], to effectively predict the signaling activity across a disease map based on the expression measures of key drug targets (KDTs). Simultaneously, it identifies the most relevant genes for accurately predicting each circuit.

Leveraging the SHAP rankings and the learner’s predictive prowess, the model implements a selection procedure to pinpoint potential drug targets (KDTs) that could regulate the disease map. This process is validated using the bootstrapping Nogueira procedure [40], which evaluates both the stability of the KDT selection method and the quality of the predictions.

The method culminates in the generation of circuit-specific confidence intervals and user-friendly thresholds for both measures, along with the identified KDTs and their associated relevance scores. For further details, please refer to Additional file 15: Methods.

Target functional prioritization

Results from the machine learning model were rescaled to a [− 1, 1] range and plotted as a heatmap using the NMF R package [41]. Calculation of the best number of KDT clusters to select was done with the Gap Statistics method (Gap Statistics chooses the number of K, where the biggest jump in within-cluster distance occurred). Then, based on the generated hierarchical clustering on the Heatmap, KDTs were divided into three groups.

Enrichment analysis of GO (Gene Ontology) Molecular Function and Biological Processes of the KDTs belonging to each of the 3 groups was done using Panther DB and GO enrichment tool [42], those terms with adjusted p-value < 0.05 were considered significant.

The circuits integrating the RP Disease Map were functionally annotated with the UniProtKB DB (accessed on 10.09.2022). From the resulting annotations, and according to the molecular pathophysiology of RP [43, 44], 9 processes were identified, grouping the functions into RP hallmarks: Apoptosis, Necrosis, Stress Response, Inflammatory response, DNA integrity, Fatty acids, and lipid metabolism, Sensory and stimuli response, Development processes, Neuronal Processes. The RP disease map circuits were arranged and functionally tagged based on these RP hallmarks, resulting in nine functional modules. Functional profiling of relevant KDTs was performed to evaluate the impact of these KDTs on the nine RP hallmarks. Since circuits are labelled with specific functions that fall under these 9 RP hallmarks, our machine learning model predicts the most relevant KDTs for each circuit. Consequently, the effects of these KDTs on the nine RP hallmarks were determined using the corresponding circuits as proxies.

Relevant KDTs were selected for experimental analyses in RP rd10 murine model based on their relevance score (belonging to the top 10th percentile), their involvement in defined functional hallmarks and in interesting pathways with RP biological and clinical relevance.

The characteristics of the drugs targeting relevant KDTs were studied according to the Anatomical Therapeutic Chemical Classification (ATC) from the World Health Organization. Fisher’s exact test was used, against the Drugbank 5.1.8 DB, to search for over-represented drug categories from general (ATC level 1 with 14 main categories) to more specific terms (ATC level 4 with 918 categories). The ATC levels 1, 2, 3, and 4 were taken into account. ATC level 5 was not considered since it is composed of the names of the compounds and does not provide any kind of grouping.

Animal model of retinitis pigmentosa

Rd10 mice were used as a model for autosomal recessive RP. These mice have a mutation in the Pde6b gene, affecting rod photoreceptor function and leading to their degeneration. In our conditions, two peaks of photoreceptor degeneration were observed at postnatal days (P)18 and P60. C57Bl/6J mice, sharing the same genetic background, served as controls. Mice were housed under specific conditions at the Research Center Principe Felipe (CIPF) and treated according to ethical guidelines (European Union Directive (2010/63/EU). Each study type (gene expression, GABA content, western blot and retinal histological quantification) involved at least eight mice per group except for histological analysis. For further details, check Additional file 15: Methods. The procedure was approved by the Committee of Ethics in Research of CIPF.

Isolation of total RNA and quantitative real-time PCR

Total RNA was isolated from frozen retinas of control and rd10 mice aged P15 to P60 (eight retinas for each group) using the NZY Total RNA Isolation Kit, following the manufacturer’s protocol. After determining RNA concentration and synthesizing cDNA, we measured the relative expression of specific genes in retinas using real-time PCR. Expression was normalized using a housekeeping gene, and then further normalized to control mice values. For further details, and a list of probes used, check Additional file 15: Methods.

Western blot

Retinas from control and rd10 mice were homogenized and resuspended in 100 μL of RIPA buffer and boiled for 5 min. Electrophoresis was carried out on 8% SDS polyacrylamide denaturing gel at 25 mA for 2 h. Samples were then transferred to a PVDF membrane and incubated with primary antibodies against GABARɑ1, KCC2 and NR1 and subsequently with secondary antibodies. Protein bands were detected using NZY Advanced ECL and quantified using the Alliance Q9 Advanced. Western blots were quantified and normalized with β-tubulin using AlphaImager 2200 (alpha innotec, Germany). Detailed information in methodology and materials used for western blot assays can be found in Additional file 15: Methods.

Retinal histology and quantification

Eyes from mice at P23 were processed for histological analysis. After fixation and cryosectioning, sections were permeabilized with 0.1% Triton X-100 for 1 h, incubated with anti-GABARɑ1, KCC2 and NR1 primary antibodies and fluorescence-conjugated secondary antibodies. After labelling and counterstaining with DAPI, the sections were mounted in Fluoromount-Gmounted. Confocal microscopy was used for observation of retinal inner nuclear layer (INL), retinal inner plexiform layer (IPL), retinal outer nuclear layer (ONL), retinal outer plexiform layer (OPL) and retinal ganglion cell layer (GCL). To confirm photoreceptor degeneration the number of nuclei at the outer nuclear layer (ONL) was measured. For further details, check Additional file 15: Methods.

Determination of GABA in retinal extracts with liquid chromatography–mass spectrometry (HPLC–MS)

GABA concentrations in retinal extracts of control and rd10 mice at P18 and P23 were measured using HPLC–MS. After homogenization, protein precipitation, and centrifugation, samples were loaded for LC–MS. Chromatography was performed, and the HPLC was coupled to a mass spectrometer. GABA concentration was determined using a standard curve and expressed in nmol/mg protein. Check Additional file 15: Methods for further details.

Statistical analysis for preclinical validation

Statistical analysis was performed using GraphPad Software 9.0 (Prism; GraphPad Software, Inc, San Diego, CA). Normal distribution of data was analyzed by Shapiro–Wilk and Kolmogorov–Smirnov tests. Comparisons between age-matched control and rd10 mice were performed using unpaired t-tests or Mann–Whitney test. A P value < 0.05 was considered statistically significant. The data were plotted using Graph Pad Software 9.0. The data were presented as mean ± SEM.

Selection of drugs targeting candidate genes

Drugs targeting experimentally validated KDTs were filtered out based on pharmacological actions that counteracts the dysregulation observed in rd10 model, which means that have the opposite effect (agonist/antagonist) on target that the observed in the RP model, both in target gene expression and protein quantification. Only drugs with potential to revert the RP dysregulation were further studied. This information was obtained from DrugBank database. Assessment of the resulting drugs over the RP hallmarks modules was carried out by measuring the percentage of the hallmark (no. of circuits from total no. of circuits of a given hallmark module) that was affected by each drug, taking the relevant KDTs that each drug targets as a proxy.

Results

Mechanistic disease map of retinitis pigmentosa

In this study, we present the mechanistic map for the study of the rare disease Retinitis Pigmentosa (RP). As a starting point, we retrieved from OrphaNet a list of 93 genes associated with RP. The gene Entrez IDs plus Gene symbols as well as the database identifiers are shown in Table 1. Since HPO is a DB that uses phenotype terms to tag genes and diseases, we have exploited this dual property to screen for the genes that share at least 10 out of 22 Retinitis Pigmentosa HPO-terms (RP-HPO), resulting in a total of 22 RP-HPO-terms, as shown in Table 2. In Additional file 11: Fig. S1A, we can see the number of genes sharing RP-HPOs added that were represented in KEGG signaling pathways. The gene symbols of the added genes per number of shared HPO-terms can be accessed in Additional file 1: Table S1 and Additional file 11: Fig. S1A.

Table 1 List of OMIM/ORPHA genes associated to RP with their corresponding Entrez ID and gene symbol

Table 2 List of the 22 HPO terms associated to RP with specificity level between 7 and 12

To maintain the specificity and not over-expand the RP mechanistic map, we have selected those genes sharing at least 10 RP-HPO-terms and that were part of signaling pathways according to KEGG DB (Table 3). Finally, the RP mechanistic map was created by merging the circuits related to the core (ORPHA/OMIM) and HPO-amplified sets of RP genes. This resulted in 226 circuits that belong to 40 KEGG signaling pathways (Table 4).

Table 3 List of genes sharing ≥ 10 out of the 22 HPO terms related to retinitis pigmentosa for extended RP-mechanistic map

Table 4 KEGG pathways involved in the Retinitis Pigmentosa mechanistic map with mapped the RP-genes

To functionally annotate the RP Mechanistic Map, we used UniprotKB annotations and the GO ontologies. From all the annotations obtained for each circuit, based on RP biological mechanisms described in the literature, we prioritized those functions related to RP progression by thorough manual curation. The GO/Uniprot functions as well as the hallmarks annotating each circuit of the RP Map are shown in Additional file 2: Table S2.

Filtered functional annotations of the RP Map converged into 9 RP related processes that we defined as the RP hallmarks: “Apoptosis”, “Necrosis”, “Stress Response”, “Inflammatory response”, “DNA integrity”, “Fatty acids and lipid metabolism”, “Sensory and stimuli response”, “Development processes” and “Neuronal processes”. A summary of the percentage of hallmark modules by pathways is depicted in Additional file 11: Fig. S1B. Pathways with a higher number of circuits, like cyclic adenosine monophosphate (cAMP) dependent signaling pathway (hsa04024), forkhead box O (FoxO) signaling pathway (hsa04068), Ras signaling pathway (hsa04014), and PI3K-Akt signaling pathway (hsa04151) show a homogeneous high coverage over the hallmarks modules. In contrast, pathways with fewer circuits (smaller in the RP map), like GnRH signaling pathway or VEGF signaling pathway, show more specificity towards specific hallmark modules like “Apoptosis” and “Development processes” (GnRH signaling pathway) or “Fatty acid and lipid metabolism” and “Inflammatory response” (VEGF signaling pathway). Likewise, the distribution of the hallmarks among the pathways is highly dependent on how broad the hallmark is. Broader hallmarks are related to many pathways compared to hallmarks that tag fewer circuits and hence fewer pathways. A summary of the percentage of hallmark modules by pathways, as well as the percentage of coverage of each hallmark on the RP mechanistic Map, are depicted in Additional file 11: Fig. S1B, C and the percentage of coverage of each hallmark on the RP mechanistic Map is depicted in Additional file 11: Fig. S1C.

The information regarding RP mechanistic map, pathways, and circuits topology is openly available at the Hipathia RP viewer, and accessible at http://hipathia.babelomics.org/RP_Mechanistic_Map/ [35].

The Hipathia RP viewer allows the user to explore the elements and filter the circuits or pathways desired. It shows the interaction among the protein nodes, the ORPHA/OMIM, and HPO amplified genes in different colors, as well as the RP-hallmark related to each circuit. The user can search for specific genes based on their gene symbol or entrez ID, or for specific hallmarks of RP from the 9 described. Figure 2 illustrates the overall layout of the Hipathia RP viewer and the RP Map.

Fig. 2

Overview of the RP mechanistic map viewer in the Hipathia Web tool. The RP mechanistic map is composed of 226 circuits belonging to 40 pathways. The Estrogen signaling pathway is shown as a selection example. Colors represent RP genes found to be associated with the disease in OMIM/ORPHANET (Blue) or to be related to the disease phenotypes in HPO db (green). The shape of the edge indicates the type of interaction (activation/inhibition) between circuit nodes. The left tab above the main panel enables the selection/search of specific proteins on the selected pathway, while the right tab allows the option to export as an SVG image the visualized pathway

Approved drug targets selection

To assess the potential effect of actionable therapeutic targets over the RP Mechanistic Map we have used DrugBank 5.1.8 DB. From 5004 targets, 7919 drugs, and a total of 26,979 KDT-drug combinations parsed, we have included 711 KDTs from 1410 approved drugs with known pharmacological action and measurement values in the GTEx V8 dataset. In total, we included 2688 combinations of drug-KDT that can be found in Additional file 3: Table S3. The drug-KDT effects were simplified into 5 types resulting in 841 drug-KDT “activators”, 1530 drug-KDT “inhibitors”, 24 drug-KDT modulators, 192 drug-KDT “ligands”, and 101 drug-KDT with various effects categorized into “other”.

Machine learning model performance

Altogether, applying the proposed drug repositioning methodology to RP, we have identified 109 KDTs, belonging to 284 drugs (Fig. 3A), whose expression patterns could potentially modulate the activity of the RP mechanistic map circuits. All the KDTs × RP circuits model scores can be found in Additional file 4: Table S4. The selection of relevant KDTs is depicted in a boolean matrix in Additional file 5: Table S5. Results from the model performance are shown in Fig. 3B and Additional file 6: Table S6 as the 95% confidence intervals for the mean R² and the Nogueira stability estimate [40] for each specific circuit. Since the trustability of the results relies on stability across the ML model predictions, 207 circuits with a stability score above 0.4 (the empirical lower effect size threshold proposed in [40]) were filtered from the 226 RP mechanistic map circuits.

Fig. 3

Retinitis pigmentosa mechanistic model from the inside. A Overview of the Retinitis Pigmentosa machine learning-driven approach to drug repurposing. A Multi-output random Forest (MORF) regressor was trained, with gene expression from the GTEx dataset, and the circuit’s activity values, to assess the predictive power of 711 KDTs over the activity of the 226 circuits that configure the RP mechanistic map. Feature selection of highly predictive KDTs for a given circuit was calculated based on the R2 score and up to the top 5% KDTs ranked by the mean absolute Shapley value. B Machine learning model (MORF) stability vs R2 score circuit metrics. C Boxplots representing the distribution of the absolute model relevances (Y-axis) across the top 30 ranked KDT (X-axis). The rank is obtained by means of the L1 norm and the values are colored by their sign (positive = red; negative = blue)

Relevance of known-drug target genes in retinitis pigmentosa mechanistic map

The model selected 109 KDTs out of 207 circuits that compose the RP map and assigned a relevance score to each KDT-circuit pair (see “Methods”). The score shows the KDT’s impact on the circuit, either positively or negatively. The rescaled score (0–1) distribution of the top 30 scoring KDTs is depicted in Fig. 3C, highlighting the positive/negative regulating role of certain KDTs which could potentially indicate which drugs, based on their effect over the KDT, would better revert the disease effects.

To exploit the full potential of the results, we generated an annotated heatmap with the relevant KDTs as columns, the stable circuits as rows, and annotations on the most frequent drug effect over each relevant KDT on top (Fig. 4A). The intensity of the color represents the strength of the KDT interaction over a specific circuit, while the sign of the interaction is shown in blue (negative) or red (positive). Furthermore, to gain insight into the KDTs scores and their structure, we used hierarchical clustering of the resulting relevance matrix and plotted the distribution of scores of each relevant KDT (Fig. 4C). We used the Gap statistic among other metrics to determine the optimal number of clusters (Additional file 12: Fig. S2A), resulting in a total of 3 clusters, cluster 1 represented by 4 KDTs, cluster 2 by 15 KDTs and cluster 3 by 90 KDTs (Additional file 7: Table S7). As we can observe in Fig. 4A, a small group of KDTs, from cluster 1, are strongly influencing the activity of most circuits. Interestingly, these 4 KDTS, are either part of GABA and glutamate receptors (GABRA1 and GRIN1), related to important neuronal-development processes, involved in chemical synaptic transmission (HRH3), or involved in ion transmembrane transport processes (SLC12A5). These processes, which play important roles in various aspects of retinal function, were also highlighted by the results obtained in the Gene Ontology (GO) enrichment analysis (Additional file 12: Fig. S2C) of cluster 1. On the other side, KDTs belonging to clusters 2 and 3 seem to influence only certain regions of the RP mechanistic Map, showing a higher specificity. Results from the GO enrichment analysis show that cluster 2 is strongly focused on ion transport and synaptic signaling processes while cluster 3, which exhibits a high variability (Additional file 12: Fig. S2B), embraces a wide range of biological processes, some of them related to tissue remodeling and cell migration. Additionally, cluster 3 is associated with processes related to cellular response to hormones and immune response, while cluster 2 is more focused on nervous system development and cellular communication. Despite the differences, all clusters share some commonalities such as processes related to the regulation of ion transport as well as to cellular response to chemical stimulus, regulation of cellular component organization, and the activity of neuronal receptors (Additional file 7: Table S7, Additional file 12: Fig. S2C).

Fig. 4

Influence and profile of KDTs on Retinitis Pigmentosa mechanistic map. A Heatmap of the 109 relevant predictive KDTs (X-axis) normalized (− 1, 1) SHAP scores over the stable circuits of the RP Map (Y-axis). The sign indicates the direction of the KDT influence over each specific circuit and the score value depicts how strong is the influence of a specific KDT for predicting the activity of a specific circuit. The top color bar represents the most frequent drug effects of that specific target (KDT). At the bottom, the KDT groups obtained by hierarchical clustering are represented by red (cluster 1 with 4 KDTs), yellow (Cluster 2 with 15 KDTs), and green (Cluster 3 with 90 KDTs) color. B) Barplot representing the distribution of SHAP scores showing the accumulated volume of positive (in red) and negative (in blue) scores (X-axis), for each specific KDT (Y-axis) overall RP Map. On the left side, KDTs are annotated with their gene symbol and colored by their cluster. C Chord diagram representing the shared influence of relevant KDTs over the circuits belonging to the 9 RP hallmarks modules. The size of the chord joining two hallmarks represents the number of shared KDTs, with their SHAP scores weighted by the % of the hallmark covered (no. of circuits tagged by the hallmark with the KDT/total number of circuits with the hallmark)

Furthermore, the distribution of KDT scores depicted in Fig. 4B highlights the dual role (positive = activating; negative = repressing) that proteins might exhibit depending on the cellular process they are involved in. Although most KDTs show a balance between positive and negative scores, others like FGF2 and ROCK1 seem to have a higher proportion of positive scores condensed in fewer circuits, which are mostly related to Development and Neuronal Hallmarks. Due to the proportions of the results, only clusters are shown in the Heatmap. An explorable version of the Heatmap with circuit names and KDT symbols can be accessed at Additional file 13: Fig. S3.

In order to assess how the RP hallmarks were being modulated by the KDTs, we identified the shared KDTs among the circuits belonging to each of the RP hallmark modules, giving weights to the interaction based on the KDT mean score. As observed in Fig. 4C, there is a strong link between the “Neuronal processes”, “Development”, and “Apoptosis” hallmark modules. Overall, since one circuit can be tagged by multiple hallmarks when aggregating the influence of KDTs, most hallmarks are influenced by a wide range of KDTs.

In our study, we utilized a functional profiling approach to select KDTs for validation in rd10 mice. By gathering the set of annotations linked to the RP map circuits into the 9 RP hallmarks, we can (functionally profile the KDTs) assess how a relevant KDT affects each functional module defined by the hallmarks. Our Disease Map, linked to nine RP-associated hallmarks, prioritized KDTs mostly affecting circuits central to neuronal and development processes. Taking into account that Retinitis Pigmentosa (RP) progression encompasses not only photoreceptor cell dysfunction but also remodeling of the inner retina and its synaptic connections, our emphasis was on neurotransmitter-mediated communication. Key to this communication are glutamate and GABA, governing vertical and lateral retinal neuron communications, respectively. Given their significance in maintaining visual system integrity and their roles in neuronal and developmental processes, we specifically selected top-scored KDTs associated with receptors of glutamatergic, GABAergic, and Glycinergic neurotransmission. These KDTs were chosen based on their prominent functional roles in these RP-related processes, providing a focused approach to our experimental validation.

Dysregulation of GABA, glycine, and glutamate receptor subunits in rd10 mice

In this study, we evaluated the overall impact of target genes for the GABA A receptor subunits alpha and epsilon (Gabra1 and Gabre), the glutamate ionotropic receptor NMDA type subunit 1 (Grin1), the glycine receptor alpha 2 subunit (Glra2), the solute carrier family 12 member 5 (Slc12a5) that encodes K+–Cl− cotransporter 2 (KCC2) in rd10 mouse retinas compared to age-matched control retinas from postnatal days P15 to P60 (Fig. 5). In addition, we previously evaluated the role of target genes lipoxygenase 5 (Alox5) and elongase 4 (Elovl4), both related to fatty acid metabolism, in the mechanisms of RP at early postnatal days [48]. We observed a significant downregulation of Slc12a5 gene (Fig. 5A, B) and a less acute downregulation of Gabra1 and Gabre genes at P15, just before the peak of rod degeneration happens to be observed (P18) (Fig. 5A, C, D). At P18, we detected an upregulation for Gabre gene accompanied by a decrease in GABA content in retinal homogenates (Fig. 5A, D, E). We detected significant upregulation of Slc12a5, Gabra1, and Gabre genes from P23 to P60 (Fig. 5A–D). The upregulation in gene expression of Slc12a5 and Gabra1 was confirmed by a western blot of their proteins KCC2 and GABARα1 (Fig. 5H, I). As shown in Fig. 6A, photoreceptor degeneration is quite evident at P23. At this age the number of nuclei remaining in the outer nuclear layer (ONL) in rd10 retinas was significantly lower than in age-matched control retinas. We also performed an immunostaining showing their localization mainly at the OPL, INL and IPL at P23 (Fig. 6B, C). Retinas from rd10 mice showed a higher content of KCC2 and GABARα1 than age-matched retinas from control mice at the inner retina at this age (Fig. 6B, C). At P18 we did not observe significant changes in the target genes Glra2 or Grin1. However, from P23 to P60 both genes were upregulated (Fig. 5A E, F. The upregulation in gene expression of Grin1 was also confirmed by western blot and immunostaining of its protein NR1 (Figs. 5J, 6D). We previously described that Elovl4 gene was dramatically decreased in retinas from rd10 mice compared to age-matched retinas from control mice from P13 to P60 (data not shown). Besides, we showed downregulation of Alox5 gene in the early stages of RP, up to P23, and from there it started to gradually increase until it was significantly upregulated at P60 (data not shown). Taking together these findings, we proposed that RP leads to an altered synthesis of specialized pro-resolving mediators which are involved in the resolution of inflammation [44] (data not shown). All data can be accessed in Additional file 8: Table S8.

Fig. 5

Experimental validation of selected KDTs in rd10 mice, murine model of Retinitis Pigmentosa. Normalized gene expression profiles of selected KDTs from postnatal day (P) 15–60 with the observed peak of rod degeneration at P18 (A). Gene expression of K+/Cl− co-transporter 2 (Slc12a5) (B), alpha 1 and epsilon subunits of GABA A receptor (Gabra1 and Gabre) (C, D), glutamate ionotropic receptor NMDA type subunit 1 (Grin1) (E), and glycine receptor alpha 2 subunit (Glra2) (F), retinal GABA concentration in retinal extracts (G, the protein content of KCC2 (H) and GABA_Rα1 (I) and NR1 (J) with representative western blots at P23 in control and rd10 mice. The unpaired t-test was used to compare control “C” and rd10 mice. *p < 0.05; **p < 0.01; ***p < 0.001. All data were presented as the mean ± standard error of the mean (SEM). We analyzed at least eight retinas for each group and age

Fig. 6

Histological study of selected targets in retinas sections of rd10 mice, murine model of Retinitis Pigmentosa. Confocal micrographs of transverse sections through rd10 mice retina. DAPI nuclear staining illustrating the retinal cytoarchitecture (retinal inner nuclear layer INL, retinal inner plexiform layer IPL, retinal outer nuclear layer ONL, retinal outer plexiform layer OPL and retinal ganglion cell layer GCL). Photoreceptor degeneration (A, B), and retinal distribution of selected KDTs KCC2 (C, D), GABA_Rα1 (E, F) and NR1 (G, H) in control and rd10 mice at P23. The Mann–Whitney test was used to compare control “C” and rd10 mice. Scale bar: 50 µm. *p < 0.05. All data were presented as the mean ± standard error of the mean (SEM). We analyzed at least four retinas for each group and age

Functional profiling of drugs targeting candidate genes

Out of 1410 drugs, 284 were found to target relevant predicted KDTs. Details of these drug-KDT interactions can be found in Additional file 9: Table S9. The distribution of drug effects is shown in Fig. 7A. Over-represented categories from the different ATC levels with an FDR-adjusted p.value < 0.5 are shown in Fig. 7B and provided in Additional file 10: Table S10. The ATC classification ranges from broad categories (ATC level 1) to specific ones (ATC level 4), excluding the compound names (ATC level 5). Figure 7B reveals an enrichment in Nervous-system and Antineoplastic-Immunomodulating drugs at ATC level 1. Deeper levels highlight a significant presence of calcium channel blocker drugs, with 86% coverage at ATC level 2 and 91% at level 3. ATC level 4 also indicates a high concentration of interferon drugs, serotonin agonists, barbiturates, and various analgesic drugs.

Fig. 7

Functional profile of potentially repurposable drugs. A Venn diagram representing the simplified pharmacological actions of drugs, collected in Drugbank 5.1.8 db, with at least one target (KDT) predicted as relevant by our model. B Over-representation analysis of drugs with relevant KDTs over the ATC categories (Y-axis) across the 4 ATC levels, showing the total number of drugs that fall into each category (X-axis) and the percentage of coverage of relevant drugs over the parsed DrugBank 5.1.8 db. Only ATC categories with an FDR p-value ≤ 0.05 are shown. Fisher’s exact test was used to identify overrepresented ATC categories in our relevant results compared to their proportions in the parsed DrugBank 5.1.8 db. C Spider plot illustrating the profile of selected promising drug candidates, Omega-3 carboxylic acid, Flumazenil, and Taurine, targeting experimentally validated targets ELOVL4, GABRA1 and GLRA2 respectively, over Retinitis Pigmentosa functional hallmarks. Drug candidates were selected based on having a pharmacological action with the potential to revert the observed RP effects in the expression patterns of their targets (relevant KDTs), measured on rd10 RP murine models

Based on experimentally validated KDTs, we identified drugs, listed in Table 5, whose effects over their targets could potentially reverse the disease effects (the UP/DOWN regulation observed in rd10 mice). These include drugs for GABA receptors, ELOVL4, and Glycine receptors. To understand the full drug impact on the RP Mechanistic Map, we measured how the drugs would influence the RP hallmark modules, based on the scores we obtained for their targets. Figure 7C shows that compounds like Omega-3 and fish oils significantly affect the “Neuronal”, “Inflammatory response”, and “Fatty Acids and Lipid processes” modules. Flumazenil, a GABA receptor antagonist, primarily impacts the “Stress Response” module. Taurine, vital for many biological functions, notably affects the “Necrosis” and “Sensori and stimuli” modules, aligning with its cytoprotective properties and high presence in excitatory tissues.

Table 5 Drugs targeting validated KDTs with the potential to revert RP effects, highlighting functionally studied promising drug repositioning candidates

Discussion

In this work, we propose a unique approach that, starting from the set of disease-affected genes, provides a comprehensive landscape of the molecular mechanisms of the RP disease along with its druggable space. The method establishes the RP mechanistic disease map as an actionable environment and employs an explainable machine-learning model to assess the influence of druggable molecules, like KDTs, over the disease environment. Our approach merges information from transcriptomics, pathway graphs, biological/clinical DBs, and drug-target interactions, to generate an in-depth view of the disease. The novelty of this workflow lies in the integration of multiple data sources, reinforcing interpretability with biological knowledge while reducing the dimensionality of the datasets.

Our methodology has enabled us to create a comprehensive yet straightforward map of RP, representing its functional landscape. The utility of this method relies on its interpretability. Artificial intelligence is often used in repurposing methods to compensate for the lack of prior knowledge. However, rational target discovery and drug repositioning can be a useful addition when prior knowledge is available, as it can help fill in the gaps and limit the unknown variables. By employing mechanistic models grounded in solid biological knowledge, we have managed to reduce dimensionality by providing a causal link between gene expression and functional cell behavior, and assess the impact of specific targets on different parts of the mechanistic map. This is particularly beneficial for diseases with diverse genotypic presentations and complex pathophysiology, such as RP [5].

Retinitis pigmentosa (RP) is a heterogeneous genetic disorder that leads to the progressive degeneration of photoreceptor cells in the retina, eventually leading to blindness. Although the pathogenesis of RP is complex, as part of this work, we have identified key features of the disease. Based on the signaling circuits identified by our methodology, and supported by the scientific literature, we have defined nine hallmarks of RP: “Apoptosis”, “Necrosis”, “Oxidative stress”, “Inflammatory response”, “DNA integrity”, “Neuronal processes”, “Fatty-acid–lipid metabolism”, and “Sensory and stimuli response”.

Photoreceptor cells, like other neuronal cells, have a postmitotic nature that allows for extended cell life beyond normal physiological conditions, but it may also induce abnormal photoreceptor cell death [45, 46]. We have indeed defined cell death (“Apoptosis” and “Necrosis”) as RP hallmarks within the proposed RP actionable map.

In RP, apoptosis is the major mechanism of photoreceptor cell death, involving the activation of caspases or certain members of the TNF receptor family [47], triggered both by intracellular stress signals, such as DNA damage or oxidative stress, processes defined as RP hallmarks, or by the binding of extracellular ligands to death receptors on the cell surface [43, 48].

Necrosis also has been linked to retinal degeneration as the main mechanism responsible for massive uncontrolled photoreceptor cell death [49], necrosis-like processes can happen rapidly and are linked to a failure in bioenergetics [46]. Interestingly, receptors that initiate cell death are part of the TNF superfamily, a family of cytokines involved in inflammation (another defined hallmark of RP) and in several pathways within our defined RP map [50, 51]. Inflammation is usually mediated by TNFα, IL-6, IL-1α, and IL-1β, among other mediators, all of which are present in RP patients and murine models of the disease [52,53,54], regulating the production of chemokines and cytokines [55,56,57,58]. Moreover, TNFα can promote, via RIPK3-MLKL signaling, the activation of NLRP3 inflammasome, a major multiprotein that induces inflammatory-mediated immune cell infiltration, accelerating necrosis-like photoreceptor cell death [58,59,60].

Oxidative stress is another critical component of the pathogenesis of RP that is present in both the RP map and the predicted relevant KDTs. In RP, oxidative stress is initiated by the release of ROS from dying photoreceptor cells and microglia cells. The accumulation of ROS leads to the activation of additional stress pathways, including the unfolded protein, leading to cell death [44, 60,61,62].

We also observed the implication of energetic and glucose-related functions, including HIF1 and HK1 pathways. Both HKI and HKII isoforms are expressed in the retina and target genes of the HIF-1 pathway, part of the RP map [58]. In fact, several cases of mutations in the HK1 gene were reported to cause RP in humans [4, 5]. This could affect glucose metabolism and energy supply to retinal cells, given that photoreceptor cells have very high energetic requirements.

Since neurons are highly excitable and sensitive to ischemia-induced death, due to their high ATP turnover rate [48, 63], the presence of voltage-gated sodium or calcium channels makes neurons susceptible to overload and swelling, while glutamate-gated channels make neurons susceptible to excitotoxicity and rupture [64,65,66]. Thus, stimuli and sensory transduction also plays a role in RP progression, making photoreceptor cells and other neural cells susceptible to excitotoxic neurotransmitters, such as glutamate, or to the action of other neurotransmitters, like GABA or glycine [67]. As mentioned above, RP leads to a progressive loss of photoreceptor cells, but other neural cells suffer morphologic and metabolic changes resulting in retinal remodeling including loss of retinal ganglion cells (RGCs) or network rewiring [68].

As commented above, gene defects primarily lead to the dysfunction and death of photoreceptor cells and many therapeutic approaches are focused on these cell types (outer retina). However, as RP progresses a remodeling of the inner retina also occurs including formation of abnormal neural circuits and changes in the electrophysiologic properties of the retinal network [69,70,71]. After photoreceptor cell death, the inner retina changes their synaptic connections affecting the transmission of information [72]. The communication between retinal neurons is mediated electrical but also chemical signals, the former includes the neurotransmitters glutamate, and GABA and their receptors. Glutamate mediates vertical communication between photoreceptor-bipolar cells and bipolar-ganglion cells. While GABA mediates lateral communication via horizontal and amacrine cells. The maintenance of the integrity of the visual systems beyond photoreceptors is a limitation for the successful outcome of any retinal treatment. In this study, we selected as KDT some subunits of the receptors associated to the glutamatergic, GABAergic and Glycinergic neurotransmission because changes in these receptors would have a significant impact in retinal cross-talk. In addition, the understanding of the time-course of the retinal remodeling may be important for achieving a proper window of therapeutic intervention [69].

In this study, we predicted and experimentally validated that RP is accompanied by alterations in neurotransmitter subunits, such as glutamate, GABA, and glycine. In many types of horizontal and amacrine cells, GABA serves as the main inhibitory neurotransmitter, in addition to glycine, which is known to also have an inhibitory role [73]. Retinal activity functions required a balance between excitation (glutamate) and inhibition (GABA and glycine). By increases in chloride conductance, GABA causes membrane hyperpolarization. Active chloride extrusion is achieved via the coupled transport of potassium and chloride ions via K–Cl cotransporters (KCC), with KCC2 being the primary active chloride extrusion system responsible for GABAergic systems. Hence, the changes in KCC2 expression (predicted and validated) could affect the GABAergic driving force [74]. Indeed, the ionotropic GABA A receptors (consisting of five subunits) are predicted as relevant KDTs, and found in almost all retinal cells [75, 76]. Inhibitory neurons seem to influence how the excitatory neurons integrate information through synaptic interactions. Under physiological conditions, retinal inhibition modulates the threshold responses of RGCs under ambient light, selectively masking certain signals to ensure the appropriate ones cross the optic nerve. GABA A receptor seems to be involved in the dynamic control of this masking inhibition and serves as a mechanism for neuronal adaptation. The upregulation of the alpha 1 and epsilon subunits of the GABA A receptor (Fig. 5C, D) observed in rd10 retinas could potentially change the threshold of RGCs activities leading to modifications in the perceived visual environment [77].

RP leads to photoreceptor degeneration triggering changes in the retinal morphology in human and animal models of RP. Studies on rd1 mice, carrying a mutation on the exon 7 of the Pde6b gene encoding the beta subunit of cGMP-PDE, with a greater rate of degeneration than rd10 mice, show that cone-mediated GABAergic amacrine cells displayed functionally altered NMDA receptors (such as Grin1) and high levels of GABA in rod bipolar cells, processes predicted by our model and validated as significantly affected in rd10 mice (Figs. 5E–J, 6 and Additional file 14: Fig. S4). In other studies, the use of GABA receptor antagonists in P23H rats (RP animal model carrying a missense mutation in the rhodopsin (rho) gene) reduces the observed stimulation thresholds of RGCs [78, 79]. Interestingly, the use of picrotoxin, an antagonist for GABA A/C and glycine receptors, also results in the blockade of the inhibitory mechanisms and ameliorates visual function in rd10 mice [80].

Although the role of GABAergic systems in RP is not completely established, our data suggest that the GABAergic system may be overactive during RP progression. We observed a significant increase in Gabra1 and Slc12a5 genes expression, their respective proteins GABARα1 and KCC2 as well as a similar distribution of them in the outer and inner plexiform layers (OPL and IPL) in retinal sections at P23. Functionally, GABA A receptor and KCC2 are both membrane proteins, closely related, which modulate GABAergic inhibition. An upregulation of Slc12a5 expression could increase Cl- extrusion, and subsequent Cl− gradient through the membrane facilitating the entry of Cl− through the ionotropic GABA A receptor. In addition, we showed an upregulation of Gabra1, which could exacerbate the inhibition induced by GABA through its binding to GABA A receptors. The higher staining of KCC2 and GABARα1 observed in the inner retina of rd10 mice (Fig. 6B, C) would support this idea. In addition, it has been suggested that the remodeling of the inner retina is accompanied by upregulation of glycine and GABA A receptors in rd1 mice [81]. Our findings suggest that therapeutic approaches aimed at reducing GABAergic inhibition and increasing the excitation of RGCs may be beneficial for RP. It has been suggested that anticonvulsants, such as tiagabine and vigabatrin, have a protective role in light-induced retinal degeneration in Abca4^−/−Rdh8^−/− mice, and that this protection is mediated by their role as GABA modulators [82]. Novel antiepileptic drugs act selectively through the GABAergic system and, although the underlying mechanism relating antiepileptic drugs with retinal degeneration is not yet fully explored, our model highlights the role of anticonvulsants in the mechanistic map of RP, probably through GABA-related processes [83]. Digging into the drugs targeting relevant predictive KDTs, the overrepresentation of anesthetic, hypnotics, and sedative drugs such as Pentazocine, Ketamine, Taurine, Flumazenil, and Phenobarbital among other drugs predicted by our model, coherently remark the relation between the anti-excitotoxic effects of certain anesthetics and their neuroprotection that has been long studied [84], whilst its role (and potential use) in RP remains undiscovered [85]. As a matter of fact, selected predicted drugs shown in Fig. 7C share a neuroprotective effect, also appearing to be related to fatty acid processes. Omega 3-related compounds, also predicted by our model, appear to have an impact on the progression of RP as discussed in Olivares-Gonzalez et al. [44].

The effectiveness of applying mechanistic models and ML methodologies to drug repurposing has been undoubtedly demonstrated in several works [17, 19, 86]. These approaches allow for the virtual testing of hundreds to thousands of scenarios with very minimal resources. While many ML-based repurposing strategies focus on the physicochemical properties of the drugs [87], they often overlook the overall impact on diseases or cells. This is where mechanistic models of disease come into play, offering insights not only into disease mechanisms but also into drug actions in specific disease contexts [88]. Indeed, supervised ML has been successfully used to boost how these tools model cellular behavior, improving the accuracy of mechanistic models in characterizing cell types [14, 49], identifying targets for disease treatments [17, 42], or predicting the effect of single or combined drugs on cells [50, 51].

The main purpose of this work was to detect actionable target genes that could have an impact on RP, providing an explorable resource of knowledge. Nevertheless, this methodology can also be used to model all cellular mechanisms to evaluate potential undesired off-target and side effects as well as to the identification of new disease-related genes, such as GABRA, GRIN1, KCC or SLC12A5 whose mutations can lead to the development of retinal dystrophies. Due to its exploratory nature, the lack of ground truth in the validation process comes as a limitation of the method, we have demonstrated the robustness of the results experimentally validating the role of selected KDTs, with no previous association with RP, in RP murine models. However, further research is needed to assess the in vitro and in vivo impact that the drugs targeting relevant KDTs would have on the disease progression or initiation.

Nevertheless, given that RP has more than 250 genes identified, but only 50% of RP cases are attributable to identified genes, this application on its own has tremendous potential for the clinical community, as a tool to evaluate the impact of genes of uncertain significance [5].

Conclusions

In this work, an actionable map of Retinitis Pigmentosa (RP) is presented, reflecting the mechanisms and dynamics of the functions involved in disease initiation and progression. The map serves as a resource for further analysis, validations, and computational modeling, such as network or graph analyses. The methodology identifies potential targets impacting RP processes, providing a list of relevances and pinpointing the mechanism, function, and circuit affected as well as the KDT responsible for the predicted impact. These relevances can guide and aid in finding possible clinical intervention points. Although focused on signaling pathways, the model can be expanded to other types of pathways using directed interaction networks, and various datasets representing gene/protein abundance, bearing in mind that limitations in sample numbers and balanced dimensionality must be considered. The model’s strength lies in its applicability to scenarios with limited knowledge, such as rare or emerging diseases like COVID-19 [19, 21, 86]. It can be adapted to infer interactions between genes and mechanisms, providing a sophisticated way to model cell behavior under disease conditions. Combining mechanistic models with ML provides a useful tool to explore disease druggable space, mechanisms, and pathogenesis aiding in genetic diagnosis and identifying genes with a key role in diseases. The possibility of using the available genomic repositories to formulate hypotheses and extrapolate conclusions in the rare diseases field constitutes an attractive opportunity to guide and accelerate the experimental validations in order to focus and prioritize target candidates.