Exploring the Molecular Targets for the Antidepressant and Antisuicidal Effects of Ketamine Enantiomers by using Network Pharmacology and Molecular Docking

1. Introduction

Major depression disorder (MDD) is the leading cause of disability and the largest contributor to suicide deaths in the world [1]. According to WHO (2017) [2], every year, about 800,000 individuals take their own lives. Chronic stress events, such as the Covid-19 pandemic, have been linked to the development of depression and anxiety as well as suicidal ideation [3]. The pathophysiology of MDD is complex and multifactorial, involving the functioning of several neurotransmission systems besides monoamines. The modulation of neuroplasticity events particularly those related to BDNF signaling, and the process of neuroinflammation are events closely implicated in the development of depressive symptoms [4,5] Despite distinct molecular mechanisms, suicide attempts (SA) and suicidal ideation demonstrate a shared genetic etiology with MDD, and severe depressive symptoms play crucial risk factor [6,7].

Ketamine is a racemic mixture of the S-ketamine (esketamine) and R-ketamine (arketamine) that has been used due to its anesthetic and analgesic properties since the 1970s. Recently, esketamine was approved by the FDA for the treatment of refractory depression and suicidal ideation because of its rapid, potent, and sustained effect, which substantially differs from the antidepressant profile of conventional antidepressants. The pharmacological properties of ketamine are mainly attributed to the inhibition of N-methyl-D-aspartate receptors (NMDAr). However, ketamine has multiple pharmacological targets, such as gamma aminobutyric acid (GABA), dopamine, serotonin, opioids, and cholinergic receptors [8]. A single subanesthetic dose of nasal esketamine produces rapid and long-lasting effects in patients resistant to other antidepressants [9]. Evidence indicates that the rapid antidepressant effect of ketamine may involve antagonism of NMDAr on GABAergic interneurons, preventing the inhibition of the glutamatergic tone, as proposed by the disinhibition hypothesis. As a consequence, AMPA receptors are activated, causing cell depolarization and activation of voltage-dependent cell-gated channels. Calcium entry then stimulates the exocytosis of synaptic vesicles containing BDNF, which causes the activation of TrKB receptors, culminating in the activation of mTOR complex 1 (mTORC1) [8,10,11].

One possibility to study the multitarget potential of ketamine is through bioinformatics tools, which allow screening of depression and suicide molecular targets rapidly and in an early approach that can guide future higher cost experimental studies. Bioinformatics is an area with great potential application for the study of diseases and drug design and discovery. Pharmacological target network analysis, molecular docking, and gene ontology (GO) enrichment analyses may allow the understanding of the polypharmacological action of molecules, investigating their potential targets and allowing a more rational use of animal models in scientific research [12,13,14,15,16]. Computational approaches are already used by the pharmaceutical industry for the development of new drugs, because besides the economic advantage, they have the potential to accelerate the discoveries [17,18,19].

In this study, we use some of these bioinformatics tools to investigate the potential molecular targets of ketamine enantiomers in depression and suicide. Ketamine has been reported to exhibit polypharmacological action, and there may be differences between esketamine and arketamine in terms of their potency and biological effects that may be due to their ability to interact with different molecular targets [20,21]. Therefore, we also computationally investigated the overlapped and/or unique molecular targets of esketamine and arketamine in depression.

2. Results

2.1. Identification of the ketamine enantiomers, suicide, and depression-overlapping genes

According to SwissTargetPrediction database, 105 potential targets of arketamine and 105 potential targets of esketamine were found. Importantly, although most of the target genes are shared by the enantiomers (91 genes), there are differential genes that are unique targets for arketamine (14 genes), and others unique to esketamine (14 genes).

As a result of the search for depression-related genes and suicide-related genes in the GeneCards database we obtained a total of 13,914 genes and 2,409 genes, respectively. Of note, suicide and depression share 1983 genes, according to the prediction.

Venn diagram illustrated in Figure 1 shows potential targets among depression-related genes, suicide-related genes, esketamine and arketamine-related genes. The intersection of arketamine and esketamine, depression and suicide-related genes yielded 53 potential targets. Venn diagram also shows 31 genes intersected by depression, arketamine and esketamine, 6 genes intersected by depression, suicide and arketamine (PDE4B, KDM4C, MIF, NOS2, CREBBP, HMOX1), and 6 genes intersected by depression, suicide and esketamine (SIRT1, PIK3CA, ICAM1, GABRA1, GABRB2, PRF1).

2.1.1. Network construction of ketamine enantiomers, suicide, and depression-overlapping genes

The search for the ketamine enantiomers, suicide and depression overlapping genes in the Homo sapiens organism was performed In GeneMANIA. The network obtained was imported and analyzed in Cytoscape V. 3.9.1. software. As a result, a network with 73 nodes, 891 connections, and a clustering coefficient equal to 0.349 was obtained. The colors of the edges (connections) reflect the type of interaction analyzed. The network has 891 connections, in which 34.85% are physical interactions (green), 28.20% are co-expression (yellow), 25.94% represent shared protein domains (pink), 9.09% are co-localizations (blue), 7.20% are predictions (gray), 2.84% are genetic interactions (purple), and 1.88% are metabolic pathways (black). The network construction of esketamine, arketamine, suicide, and depression-operlapping genes is shown in Figure 2A. Regarding the network organized in Cytoscape software, the larger the size of the node, the higher the degree (n° of connections) of a given protein, and the darker the color, the higher the betweenness centrality (importance of that protein as a "bridge" or shortest pathway in the network). When we built the network with 53 overlapping genes, 20 genes were retrieved in the network, which are: KCNN4, NQO1, AOC2, VEGFD, VEGFC, CA3, AOC1, CA12, GABRA3, PDE5A, MTNR1B, GABRB2, CA13, GRIN2C, CA8, GABRA6, GABRA4, GRIN3A, PDE4B and PDE3B. Subsequently, the network was filtered considering degree>30 and betweenness centrality between 0.017 and 0.044. From these data, 11 genes were found with the highest number of connections (Degree>30) and with the shortest paths (betweenness centrality) in the network. The genes found were: GRIN2B, GRIN2A, GRM4, EGFR, KDR, GABRB3, GABRB2, GABRA3, KCNN1, PDE5A, and MTNR1B. The proteins encoded for these genes are: NMDA receptor (GRIN2B and GRIN2A), metabotropic glutamate receptor (GRM4), GABA receptor subunits β2, β3 and α3 subunits (GABRB3, GABRB2, GABRA3), epithelial growth factor receptor (EGFR), melatonin receptor 1β (MTNR1B), kinase insertion domain receptor (KDR), member 1A of the N subfamily of calcium-potassium activated channel (KCNN1), and the enzyme phosphodiesterase 5A (PDE5). The filtered network of targets obtained is illustrated in Figure 2B.

2.1.2. Network construction of esketamine, suicide and depression-overlapping genes

In the second step of the network analysis, we obtained potential targets shared by esketamine, suicide and depression. Therefore, the green diagram was built by searching for the genes SIRT1, PIK3CA, I-CAM1, GABRA1, GABRB2 and PRF1. As a result, we obtained a network with 26 nodes, 242 connections, and a clustering coefficient equal to 0.714. Figure 3A shows a high degree of clustering among the genes encoding GABAergic receptor subunits, which are more connected to each other than to other genes. Most interactions are physical interaction (48.92%), followed by predicted interactions (37.4%), co-expression (10.2%), and shared protein domains (3.48%). In the network construction, we retrieved 20 genes. After applying degree (DC>2) and betweenness centrality (BC) filters between 0.012 - 0.0667, we obtained 11 most important genes in the network, as illustrated in Figure 3B. Most of these genes encode GABAergic receptor subunits (GABRA1, GABRB1, GABRA4, GABRG2, GABRB2), but genes encoding sirtuin 1 (SIRT1), an integrin alpha L subunit (ITGAL), phosphatase (FOSPI3K) and kinase (PI3K) subunits, thrombopoietin receptor (MPL), transcription factor (BCL11A) and perforin 1 (PRF1) were also identified.

2.1.3. Network construction of arketamine, suicide and depression-overlapping genes

This third network represented in Figure 4A, whose nodes are in blue, represents the common targets between arketamine, depression and suicide and was assembled based on the search for the genes CREBBP, NOS2, HMOX1, KDM4C, PDE4B and MIF. As a result, we obtained a network with 26 nodes, 189 connections, and a clustering coefficient equal to 0.293. When we performed the pharmacological network analysis of only esketamine, the result was a cluster of GABA_A receptors (with a clustering coefficient equal to 0.714, where the closer to 1, the less random the nodes of a network are involved), highlighting the genes coding for some specific subunits, such as GABRA1, GABRA4, GABRB1, GABR2 and GABRG2. Most interactions are physical interaction (80.30%), followed by predicted interactions (9.48%), genetic interactions (4,28%), co-expression (3.72%), metabolic pathways (1.12%), co-localization (0.74%) and shared protein domains (0.37%). In the network construction, we retrieved 20 genes. After applying degree (DC>16) and betweenness centrality (BC) filters between 0.039 - 0.157, we obtained the 10 most important genes of the network, as illustrated in Figure 4B. Most of these genes code for signaling pathway-related proteins (CREB1, CREBBP and COPS5), enzyme-related genes (PDE4B, PDE4D, NOS2 and HMOX1), genes related to immune responses (MIF and CD74), and cytoskeleton (ACTN4).

2.1.4. Network construction of ketamine enantiomers and depression-overlapping genes

The purple network represented in Figure 5A was built by searching for the genes ROCK1 , AKR1B1, JAK2, DHODH, TYK2, PDE4C, TGFBR1, TNKS2, GABRA4, CCND1, PDE4D, OPRD1, RPS6KA1, F13A1, MTNR1B, PDE5A, IGF1R, BRD4, PDE3B, DYRK1A, JAK1, ADORA1, SYK, CSNK1G1, ADORA2B, TNKS, AURKB, EPHB3, CA13, GRK6, ADORA3, which are the 31 target genes shared by esketamine, arketamine and depression. As a result, a network with 51 nodes, 698 connections, and a clustering coefficient equal to 0.584 was obtained. This high degree of clustering is composed of several types of phosphodiesterases. Most interactions in this network are shared protein domains (42.12%), followed by co-expression (26.24%), metabolic pathways (17.70%), physical interactions (6.43%), co-localization (5.77%), genetic interactions (0.98%) and predicted interactions (0.76%). In the network construction, we retrieved 20 genes. After applying degree (DC>11) and betweenness centrality (BC) filters between 0.038 - 0.079, we obtained the 8 most important genes as illustrated in Figure 5B. Most of these genes encode phosphodiesterases (PDE5A, PDE7A, PDE3A), opioid receptor delta 1 (OPRD1), adenosine A1 receptor (ADORA1), Janus Kinase 1 (JAK1) involved in inflammation and epithelial remodeling, a subunit A of coagulation factor XIII (F13A1) with catalytic function, and a ribosomal protein S6 Kinase A1 (RPS6KA1) related to cell growth and differentiation.

2.2. GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of ketamine enantiomers, suicide, and depression-overlapping genes

We used GO and KEGG for a global evaluation of ketamine enantiomers, suicide, and depression-overlapping-genes. According to biological process (BP) analysis of GO, the top five terms were chemical synaptic transmission, anterograde trans-synaptic signaling, trans-synaptic signaling, synaptic signaling, and regulation of membrane potential (Figure 6A). Furthermore, based on fold enrichment, the one-carbon metabolic process, the GABAergic synaptic transmission, and GABA signaling pathway are also highlighted in BP (Figure 6A). The result of molecular function (MF) describing the top five terms were transmitter-gated ion channel activity, extracellular ligand-gated ion channel activity, carbonate dehydratase activity, gated channel activity, and ligand-gated ion channel activity (Figure 6B). Moreover, based on fold enrichment, NMDA glutamate receptor activity and GABA-gated chloride ion channel activity are also highlighted in MF (Figure 6B). The top five cellular components (CC) included synapse, neuron projection, somatodendritic compartment, synaptic membrane, and receptor complex (Figure 6C). According to CC, it is also possible to verify the importance of these genes as components of GABA-A receptor complex, GABA receptor complex and NMDA selective glutamate receptor complex (Figure 6C). The top five signaling pathways according to KEGG enrichment analysis were morphine addiction, nicotine addiction, neuroactive ligand-receptor interaction, Rap1 signaling pathway, and cAMP signaling pathway (Figure 6D). Furthermore, cocaine addiction, amphetamine addiction, GABAergic synapse, glutamatergic synapse, and dopaminergic synapse are also involved in the function of these genes (Figure 6D).

2.3. GO and KEGG enrichment analysis of esketamine, suicide and depression-overlapping genes

We use GO and KEGG enrichment analysis to investigate the functions and biochemical pathways of esketamine, suicide, and depression-overlapping genes. According to BP analysis of GO, the top five terms were GABA signaling pathway, chloride transmembrane transport, inorganic anion transmembrane transport, GABAergic synaptic transmission, and regulation of postsynaptic membrane potential (Figure 7A). Moreover, based on fold enrichment, cellular response to histamine, response to histamine and inhibitory synapse assembly are highlighted in BP (Figure 7A). The top five terms of MF were GABA-gated chloride ion channel activity, ligand-gated anion channel activity, GABA-A receptor activity, GABA receptor activity and inhibitory extracellular ligand-gated ion channel activity (Figure 7B). Furthermore, benzodiazepine receptor activity demonstrated a high fold enrichment in MF (Figure 7B). The top five terms of CC were GABA-A receptor complex, GABA receptor complex, chloride channel complex, GABAergic synapse, and dendrite membrane (Figure 7C). According to KEGG enrichment analysis, the top 5 pathways were nicotine addiction, GABAergic synapse, morphine addiction, retrograde endocannabinoid signaling, and neuroactive ligand-receptor interaction (Figure 7D).

2.4. GO and KEGG enrichment analysis of arketamine, suicide, and depression-overlapping genes

To investigate the biological functions and the biochemical pathways involved in the action of arketamine in depression and suicide, GO and KEGG analyses were performed. According to BP, the top five terms were positive regulation of cytokine production, cytokine production, regulation of cytokine production, positive regulation of gene expression and positive regulation of multicellular organism processes (Figure 8A). Additionally, based on fold enrichment, negative regulation of relaxation of muscle, regulation of relaxation of cardiac muscle, adenylate cyclase-activating adrenergic receptor signaling pathway involved in heart process, cAMP catabolic process and positive regulation of killing of cells of another organism were highlighted in BP (Figure 8A). The MF top five terms were 3’,5’-cyclic-AMP phosphodiesterase activity, phosphoric diester hydrolase activity, cyclic-nucleotide phosphodiesterase activity, 3’,5’-cyclic-nucleotide phosphodiesterase activity and cAMP binding (Figure 8B). Furthermore, heme oxygenase activity, nitric-oxide synthase activity and tetrahydrobiopterin binding demonstrated a high fold enrichment in MF (Figure 8B). The only significant cellular localizations in CC were chromatin, perinuclear region of cytoplasm, chromosome, voltage-gated calcium channel complex, and calcium channel complex (Figure 8C). The top 5 pathways indicated by KEGG were tuberculosis, hepatitis B, Kaposi sarcoma-associated herpesvirus infection, cAMP signaling pathway, and pathways in cancer (Figure 8D). These pathways are related to genes and proteins important in the immune response against infections and cancer.

2.5. GO and KEGG enrichment analysis of ketamine enantiomers and depression-overlapping genes

We used GO and KEGG for a global evaluation of ketamine enantiomers and depression-overlapping genes. The BP top five terms were cyclic nucleotide catabolic process, ribonucleotide catabolic process, purine-containing compound catabolic process, purine ribonucleotide catabolic process and cyclic nucleotide metabolic process (Figure 9A). Moreover, based on fold enrichment, cAMP catabolic process and cGMP catabolic process are also highlighted in BP (Figure 9A). The MF top five terms were 3’,5’ -cyclic-nucleotide phosphodiesterase activity, cyclic-nucleotide phosphodiesterase activity, 3’,5’-cyclic-AMP phosphodiesterase activity, phosphoric diester hydrolase activity and 3’,5’-cyclic-GMP phosphodiesterase activity (Figure 9B). The top five CC terms included somatodendritic compartment, cell projection membrane, leading edge membrane, neuron projection membrane and neuronal cell body (Figure 9C). The top five pathways according to KEGG enrichment analysis were purine metabolism, morphine addiction, metabolic pathways, cAMP signaling pathway, and cGMP-PKG signaling pathway (Figure 9D).

2.6. Molecular Docking

To perform docking, the targets of each filtered network were selected by the parameters DC and BC plus the top five highest values of closeness centrality (CC). After selecting the top 5 biological targets from each network, molecular docking was performed according to the structures and parameters described in Table 1, which also shows the affinity values in kcal/mol. The results of the docking analyses reveal that depending on the molecular target, there may be differences in affinities between the ketamine enantiomers (Table 1). We selected as priority interaction targets those whose affinity values lie within the range of -8.102 to -7.426 kcal/mol, because they are values found for the GluN2B and GluN2A subunits of the NMDA receptor (encoded by GRIN2B and GRIN2A genes, respectively) which are consolidated pharmacological targets for both ketamine enantiomers, being considered our positive control in this docking study. Therefore, the main pharmacological targets of esketamine and the respective affinity values were: GluN2B (-8.102 kcal/mol), PDE7A (-7.811 kcal/mol), PDE5A (-7.629 kcal/mol), JAK1 (-7.609 kcal/mol), PDE3A (-7.505 kcal/mol), GluN2A (-7.461 kcal/mol) and GABRG2 (-7.448 kcal/mol). The main pharmacological targets of arketamine and the respective affinity values were: PDE3A (-8.120 kcal/mol), MIF (-8.008 kcal/mol), GluN2B (-7.974 kcal/mol), PDE7A (-7.826 kcal/mol), GABRB2 (-7.740 kcal/mol) PDE5A (-7.740 kcal/mol), JAK1 (-7.563 kcal/mol) and GluN2A (-7.426 kcal/mol).

Figure 10 illustrates three-dimensional images from molecular docking simulations, indicating the interaction between the target proteins and the ketamine enantiomers at the respective binding sites and the presence or absence of hydrogen bonds.

In the computational prediction, MIF (1GCZ structure) has one of the highest docking scores for arketamine (-8.008 kcal/mol), and interacts in the center of the protein, as illustrated in Figure 10A.

According to docking scores, arketamine is predicted to have a higher affinity for PDE3A and PDE5A (-8.120 kcal/mol and -7.740 kcal/mol, respectively) than esketamine (-7.505 kcal/mol and -7.629 kcal/mol, respectively), while the affinity of esketamine and arketamine is similar for PDE7A (-7.811 kcal/mol and -7.826 kcal/mol, respectively). In molecular docking simulations, we found hydrogen bonds interactions between PDE3A Leu910 and esketamine (Figure 10B), PDE7A Asp362 and arketamine (Figure 10D), PDE7A His212 and arketamine (Figure 10E), PDE5A His613 and esketamine (Figure 10F), and PDE5A His613 and arketamine (Figure 10G). No hydrogen bridges were found between PDE3A and arketamine (Figure 10C).

We also found that the GABRB2, GABRB3 and GABRA3 genes are potential targets shared by both enantiomers of ketamine. In molecular docking analyses, the proteins obtained from PDB were 6D6T (GABA_Aα1β2γ2 receptor), 6DW1 (GABA_Aα1β1γ2 receptor) and 7QNA (GABA_Aα4β3γ2 receptor).

According to computational prediction, ketamine enantiomers have a higher affinity for the 6D6T structure than 6DW1 and 7QNA structures. In addition, the chloride channel (-7.747 kcal/mol for arketamine), the benzodiazepine site (-7.448 kcal /mol for esketamine) and the neurosteroid site (-7.371 kcal/mol for esketamine and -7.387 kcal/mol for arketamine) are potential targets for interaction with ketamine. Figure 10H illustrates the 6D6T structure, which is a GABA_Aα1β2γ2 receptor, showing the neurosteroid (in left) and benzodiazepine (in right) binding sites. In molecular docking simulations, we found two hydrogen bonds interactions between GABA_A receptor Arg141 and esketamine in neurosteroid binding site (Figure 10I), and a hydrogen bond between GABA_A receptor Asn60 and esketamine in benzodiazepine binding site (Figure 10J).

3. Discussion

The main purpose of this study is to explore the potential molecular targets of ketamine enantiomers in suicide and MDD by employing network pharmacology and molecular docking simulations. Our workflow involved a sequential analysis composed by the following steps: a) identification of molecular targets of ketamine enantiomers in SwissTargetPrediction; b) collection of genes targets of MDD and suicide in GeneCards; c) creation of a Venn diagram to obtain the overlapping genes; d) construction of the network of overlapped genes; e) implementation of network filter parameters to obtain the highest closeness, betweenness and degree centrality; f) running the molecular docking simulations. Our bioinformatics analysis uncovered that both enantiomers of ketamine, esketamine and arketamine, share biological targets that may be implicated in their antidepressant and antisuicide properties. However, due to the large number of pharmacological targets found in the network analyses, we limited the molecular docking analyses to the nodes that had the best topological parameters. The main targets obtained are NMDA receptors (GluN2 subunits), GABA receptors (GABRA3, GABRB2, GABRB3 subunits), macrophage migration inhibitory factor (MIF) and various subtypes of phosphodiesterases (PDE3A, PDE7A, PDE5A). However, it is important to consider that other targets not analyzed may be relevant to the antidepressant and antisuicidal effects of ketamine enantiomers.

GO and KEGG functional analyses pointed out relevant targets for both ketamine enantiomers and were useful to direct the subsequent studies. In GO enrichment analysis, related to molecular function, NMDA receptor activity was pointed out as one of the main targets. In agreement, ketamine is a well-known noncompetitive open-channel NMDAr antagonist [8]. The GRIN gene family encodes three classes of NMDA receptor (NMDAR) subunits: the glycine-binding GluN1 (product of GRIN1), glutamate-binding GluN2 (GRIN2A, GRIN2B, GRIN2C, and GRIN2D), and the glycine-binding GluN3 subunit (GRIN3A and GRIN3B). Regarding GluN2, four different subunits (GluN2A-D) are present throughout the central nervous system [22]. Of note, NMDA receptors are composed of two obligatory GluN1 subunits, and a combination of two additional GluN subunits (GluN2A-D, GluN3A,B) [23].

We found by docking analyses that GluN2 subunits of NMDAr are potential targets for the action of both ketamine enantiomers, particularly GluN2B. Interestingly, a study by Zhang and co-workers (2021) [9] reported that both enantiomers of ketamine interact at NMDAr in the transmembrane portion of the channel at key amino acid residues, such as Leu642 in GLuN2A, Leu643 in GluN2B and Asn616 in GluN1. These amino acid residues form hydrophobic interactions and hydrogen bonds with ketamine, and mutations in these residues can decrease the potency of ketamine in blocking NMDAr. Interestingly, several studies have indicated that the antagonism of ketamine at NMDAr is implicated in its rapid antidepressant response, particularly at NR2B subunit [24,25]. The genetic deletion of GluN2B from cortical neurons caused a marked antidepressant-like effect in mice. In addition, ketamine was devoid of antidepressant-like effect in the tail suspension test in the animals that had such GluN2B deletion [24]. Further reinforcing the role of GluN2B for the antidepressant effect of ketamine, a study by Tang et al. (2020) [26] showed that CaMKIIα affects phosphorylation and extrasynaptic localization of GluN2B. Indeed, ketamine was proposed to selectively block extra-synaptic GluN2B-containing NMDARs, causing a subsequent activation of mTORC1, which in turn induces protein synthesis [27]. The docking analyses indicated a slightly higher affinity of esketamine enantiomer at GluN2B than arketamine.

Although NMDAr is the key target of esketamine and arketamine, other targets may be implicated in their antidepressant property. Indeed, it has been discussed whether the inhibition of NMDAR is critical for the mood enhancing effects of ketamine enantiomers, and a hypothesis that has been raised is that their antidepressant and antisuicidal effect may be dependent on their interaction with multiple molecular targets [28]. Regarding this issue, our study may contribute to identify potential protein targets.

Network construction of arketamine, suicide, and depression-overlapping genes identified macrophage migration inhibitory factor (MIF) as a central node. Subsequently, our study identified MIF, to our knowledge for the first time, as a potential molecular target of arketamine by molecular docking. MIF is a cytokine originally characterized by its ability to prevent migration of macrophages that is now recognized as a key mediator of the innate immune system and adaptive immunity [29]. It is a glycoprotein expressed in many different types of cells and tissues, including brain regions such as hippocampus, and cerebral cortex (regions implicated in mood regulation) [30,31]. The release of MIF, by either anterior hypophysis or monocytes/macrophages, is stimulated by inflammatory stimuli and/or stress [29]. More recently, evidence has pointed out that MIF is involved in assembly and activation of the NLRP3 inflammasome [32]. The role of MIF in MDD has been increasingly investigated in recent years. The genetic deletion of MIF caused anxiety- and depression-like behavior and impairments in hippocampal dependent memory in mice [30]. Interestingly, the blood levels of MIF were shown to be increased in patients with MDD [33]. In a study carried out with healthy university undergraduates, MIF levels were 40% higher in the individuals that presented high-depressive symptoms as compared to those with low depressive symptoms [34]. In addition, MIF levels were higher in depressive patients who later failed to respond to antidepressants as compared with those responsive to treatments [35]. Despite these data implicating MIF in the pathophysiology of MDD, the effects of ketamine on MIF levels both in clinical studies and in animal models of depression have not been explored, deserving future studies.

Another interesting finding of the bioinformatic analysis is that 3',5'-cyclic nucleotide phosphodiesterases (PDE) may be molecular targets of ketamine enantiomers. PDEs are classified in 11 subtypes based on sequence homology, structural and enzymatic features mainly related to sensitivity to modulators [36]. These enzymes are responsible for the hydrolysis of the intracellular second messengers cyclic adenosine monophosphate (cAMP) and/or cyclic guanosine monophosphate (cGMP), which are essential for synaptic plasticity, learning and memory [37]. In addition, they play an important role in mood regulation [38]. Network pharmacology indicated the possible implication of several PDE genes as targets for ketamine enantiomers in MDD, as indicated by the high degree of clustering. In spite of the fact that other PDEs have been indicated by network pharmacology approach, we selected four isoforms considering the highest centrality parameters obtained in the network analysis. Subsequently, based on the binding energy of the protein-ligand complex, PDE3A, PDE7A and PDE5A were chosen for the tridimensional visualization and analysis of hydrogen bonds. Particularly, both esketamine and arketamine have high affinity for these PDE isoforms.

PDE3 occurs as two groups of isoenzymes (PDE3A and 3B) that are involved in the hydrolysis of cGMP, and mainly cAMP [39,40]. Although role of PDE3A in the pathophysiology of MDD is not well established, this enzyme may be implicated in depressive symptoms, since cilostazol, a PDE 3A inhibitor, has been reported to decrease the severity of depressive symptoms assessed by Hamilton depression rating scale scores in poststroke depression patients [41].

PDE7 comprises 2 subclass members: PDE7A and B. This enzyme, abundant in the brain, is responsible for the hydrolysis of cAMP. The highest levels of PDE7A have been reported in the olfactory bulb, olfactory tubercle, hippocampus, cerebellum, medial habenula nucleus, pineal gland, area postrema, and choroid plexus [42]. Of note, PDE7A inhibitors have been shown to exhibit procognitive and antidepressant properties [43].

PDE5 hydrolyzes cGMP, which is a critical component of the cGMP-protein kinase G (PKG) intracellular signaling in neurons There are three isoforms of PDE5A, encoded by a unique gene, namely PDE5A1, PDE5A2, and PDE5A3 [44]. It has been reported that PDE5A single nucleotide polymorphism (SNP) may be associated with MDD [45]. The involvement of PDE5A in MDD is reinforced by the antidepressant-like effect elicited by the PDE5A inhibitor sildenafil in chronic unpredictable mild stress [46].

Considering the involvement of the above mentioned PDE in the pathophysiology of MDD, the possible role of ketamine in the modulation of the activity of these enzymes is poorly explored and needs to be further investigated.

Network pharmacology, GO, KEGG analysis and molecular docking pointed out that GABA genes are also potential targets for both enantiomers but mainly for esketamine. The GABA receptor (GABAr) family is divided into metabotropic or ionotropic receptors. GABA_Br are the heterodimeric metabotropic receptors and are associated with potassium and calcium channels via G protein. GABA_B receptors are members of class C G protein-coupled receptor (GPCR) family. These receptors have R1 subunits where orthostatic ligands bind (encoded by GABBR1 gene) and R2 coupled to G protein (GABBR2) [47]. GABA_Ar are pentamer chloride channels and have 19 genes that encode for subunits of the GABA_A receptor, which is organized into classes (α, β, γ, ρ, θ, ε, π and δ) and isoforms (α1-6, β1-3, γ1-3 and ρ1-3) [48].Thus, the GABRA1 gene encodes the GABA Type A Receptor Subunit Alpha1(GABA_Aα1) protein, the GABRB3 gene has GABA Type A Receptor Subunit Beta3 (GABA_Aβ3) as its product, and the GABRG2 gene encodes the GABA Type A Receptor Subunit Gamma2 (GABA_AG2) protein. The assembly of GABA_Ar as hete-ropentamers produces complex heterogeneity in their structure, and different subunit composition is possible. A typical GABA_A receptors are made up of two alpha subunits, two beta subunits, and one gamma subunit, and the interaction sites lie at the connections between these subunits. The GABA interaction site is located between the interfaces of the alpha and beta subunits, the neurosteroid site is located in the transmembrane domain between the alpha and beta subunits, and the benzodiazepine site is located at the interface between the alpha and gamma subunits [49]. In our study, we found that the GABRB2, GABRB3 and GABRA3 genes are potential targets shared by both enantiomers of ketamine, while GABA_AG2 coded by GABRG2 is predicted to mainly interact with esketamine.

The benzodiazepine site of GABAr is closely implicated in anhedonic behavior. Anhedonia is one of the symptoms of major depressive disorder that is closely related to suicide risk [50]. Of note, ketamine has been reported to attenuate symptoms of anhedonia, and this property might be implicated in its ability to prevent suicidal ideation [51]. A study investigating MRK-016 (a negative allosteric modulator of GABA_Ar) reported that this compound recovered hedonic behavior in mice, and that this behavior was prevented by the use of the benzodiazepine antagonist flumazenil, suggesting the involvement of the benzodiazepine site of GABA_Ar in anhedonia symptoms [52]. A clinical study showed that esketamine improved anhedonic symptoms in patients with unipolar and bipolar depression [53]. These data support the computational prediction that the benzodiazepine site of GABA_Ar may be involved in the mechanism of action of ketamine in improving anhedonic symptoms in MDD.

Neurosteroid analogs, molecules that bind to the allosteric site of neurosteroids on GABA_Ar, have potential use as antidepressants and anesthetics [54]. Reinforcing the role of GABA_A receptor to rapid antidepressant responses, brexanolone and zuranolone, positive allosteric modulators of neurosteroid site of this receptor elicit rapid and sustained antidepressant effects in postpartum depression and MDD [55,56,57].

Ketamine is known for its dissociative properties and abuse potential [58], used even as an animal model of schizophrenia at doses higher than the ones that elicit antidepressant effect [59]. In agreement, the KEGG and GO analyzes indicated that both enantiomers have the potential to activate addiction-related pathways such as nicotine, morphine, amphetamine, and cocaine addiction, due to the involvement of GABAergic and glutamatergic receptors. On the other hand, ketamine has the potential to treat many addictions [60], requiring further research in this field.

Although our study suggests novel and poorly explored proteins as possible targets for the antidepressant and antisuicidal effects of ketamine enantiomers, future experimental studies are needed to confirm the role of these targets for the effects of ketamine enantiomers.

4. Materials and Methods

4.1. Identification of the potential target genes of ketamine enantiomers in suicide and depression

The list of biological targets of arketamine and esketamine was obtained from SwissTargetPrediction [61] by searching the structure of ketamine (Pubchem CID 3821) [62], with selection of the enantiomers using the "calculate stereo" tool by Marvin JS. The lists of genes related to depression and suicide were obtained from the GeneCards database [63] by searching for the words "depression" and "suicide", respectively. With Venny v. 2.1.0, we obtained the overlapping target genes between the groups.

4.2. Gene-gene Interaction Network Construction and Analysis

The list of genes obtained from the crossover was searched in GeneMANIA [64] and a network was obtained. The search was performed only in Homo sapiens organism. The network was imported and analyzed in Cytoscape V. 3.9.1 [65]. Topological network analysis was performed to identify the most important nodes (genes) using degree, betweenness centrality and closeness centrality criteria, whose cut-off value was defined according to each network's characteristics.

4.3. GO and KEGG enrichment analysis

GO and KEGG enrichment analysis was performed in ShinyGO v.0.76.3 [66] for each network built in GeneMania and Cytoscape. In this analysis, the terms biological process (BP), molecular function (MF), and cellular composition (CC) were used for GO and KEGG, and using the organism Homo sapiens, false discovery rate (FDR) <0.05 and 20 most important pathways as filters. The lowest FDR values were used to select the most important KEGG biochemical pathways or GO analyses.

4.4. Molecular Docking

The molecular docking analysis was carried out to predict the binding affinity between arketamine or esketamine and the potential biological targets using the DockThor 2.0 software [67], in a similar way as described previously [68]. The protein database (RCSB — Protein Data Bank) [69] was used to obtain the structures of biological targets and the ligands structures of esketamine (JC9) and arketamine (RKE). All binders underwent geometric optimization in the Chimera 1.14 software [70]. The preparation of target proteins included the removal of water molecules, the addition of hydrogens, and the assignment of Gasteiger charges. In some structures it was necessary to delete the ligands/antagonists/agonists present at the binding site as well as chains or parts of chains that were not part of the grid box, to optimize the docking analysis. The PDB structures used, and the binding sites are described in Table 1. Additionally, the poses of ligands were also observed using the Chimera 1.14.Interventionary studies involving animals or humans, and other studies that require ethical approval, must list the authority that provided approval and the corresponding ethical approval code.

5. Conclusions

In summary, our study uncovered that ketamine had a polypharmacological profile with many biological targets involved in suicide and MDD. The findings indicated GRIN2A, GRIN2B, PDE7A, PDE5A and PDE3A as the most likely targets predicted for both ketamine enantiomers. Regarding S-ketamine, it was predicted that the main targets are GABRB2 and GABRG2. In the case of arketamine, it was predicted that MIF is the main target. These new targets pointed out in our study have not yet been explored with regard to their involvement in the antidepressant and antisuicidal effect of ketamine, and therefore deserve attention in future experimental studies.