A Multifaceted Computational Approach to Identify PAD4 Inhibitors for the Treatment of Rheumatoid Arthritis (RA)

1. Introduction

Citrullination (or deamination) converts peptidyl-arginine side chains to peptidyl-citrulline. Peptidyl arginine deiminase (PAD) enzymes catalyze this modification [1,2,3]. PAD enzymes serve as post-transcriptional modifiers, producing citrullinated amino acids from the arginine amino acid chain. The reaction is carried out mostly in damaged & inflamed tissues [4,5]. Peptidyl arginine deiminase 4 (PAD4), one PAD isozyme, contributes to multiple cellular functions, including controlling gene expression, programmed cell death, immunity, and cell division & differentiation, neutrophil extracellular traps production, and preventing tumor development while any disruption in its function or off-site activity contributes to the development of numerous diseases [5,6].

Histones serve as the structural framework for DNA in eukaryotes. Histones undergo many post-translational modifications, such as citrullination, phosphorylation, methylation, acetylation, and ubiquitination [7,8,9], and PAD enzymes play a part in these modifications.

The PAD enzymes in homo sapiens consist of five isozymes (PADs 1-4 and 6) which are activated by calcium with 50% sequence similarity between them [2,10]. PAD1 is more predominant in the epidermis and uterus. Skeletal muscle, brain, inflammatory cells, several cancer cell lines, and secretory glands are within PAD2 activity, while hair follicles and keratinocytes with PAD3. Granulocytes NETs formation and different types of cancers are dependent on PAD4. Finally, PAD6 activity is mostly in oocytes and embryos [1,2,6]. All isozymes are located intracellularly, but only PAD IV decisively participates in histone deimination. However, newer studies imply that PAD II also could deiminate histones [1,2,11,12]. PAD IV and PAD II, besides cytoplasm, are present in mitochondria and the nucleus. Other than histones, PAD enzymes can also target fibrinogen, filaggrin, and actin proteins. Citrullination of these proteins contributes to rheumatoid arthritis pathogenesis [13,14]. Neutrophil cells lysis during programmed cell death or by an immune attack on them in inflammations or infections intracellular proteins form the extracellular traps (NETs) to counter the foreign body. PAD 4 participates in this process and is then released into the extracellular fluid with the lysed cell components [15].

Patients with abnormal function of PADs tend to form autoantibodies against the abnormal citrullinated proteins. Anti-PAD4 antibodies, which are highly distinct in RA and linked to the presence of Anti Citrullinated protein antibodies (ACPA) [16], are detected in RA patients’ serum. The isoforms PAD2 and PAD4 are the most highly associated with rheumatoid arthritis (RA) and autoimmune diseases. PAD4 is essential in the etiology of cardiovascular disorders, autoimmunity, multiple sclerosis, lupus, Parkinson's disease, malignancies, Alzheimer's disease and more. Giving us a target to be considered in various diseases’ management [1,2,3,17]. Anti-PAD4 antibodies can suppress or augment the enzymes’ activity, reducing or enhancing the inflammatory burden [18,19] and destruction of nearby organs – Most reports of Anti-PAD4 antibodies are in RA.

Different Small-molecule inhibitors of PAD were studied in the literature, but none were approved for any type of disease. Examples include the irreversible pan-PAD inhibitors halo-amidines and the reversible specific PAD4 inhibitors GSK199 and GSK 484 [20]. Amidine derivative structures are like benzoyl arginine amide, a PAD substrate. Examples include Cl-amidine, o-Cl-amidine, and BB-Cl-amidine, with increasing potency [21,22,23]. GSK compounds bind to calcium-free enzymes with higher potency than the amidines with additional selectivity to PAD 4 [23,24,25,26,27].

Multiple studies conducted on several inhibitors faced different problems regarding safety and pharmacokinetics, while the inhibitory activity is retained in vitro and in vivo.

Table 1. Previous investigational compounds potencies.

Irreversible amidines (Pan inhibitors)	IC50 on PAD4
F-amidine	1.9µM
Cl-amidine	22 µM
Reversible compounds (Selective)
GSK199	200 nM
GSK484	50 nM

Table 1. Previous investigational compounds potencies.

Irreversible amidines (Pan inhibitors)	IC50 on PAD4
F-amidine	1.9µM
Cl-amidine	22 µM
Reversible compounds (Selective)
GSK199	200 nM
GSK484	50 nM

Few used medications have a little inhibitory activity on PAD 4 like Streptomycin, chlortetracycline, minocycline, and Paclitaxel while the investigational chloramidines are more potent than them [28,29,30]. Another problem with these medications is their known high toxicity which makes them unfavorable to further study.

Figure 1. Known substances with PAD 4 activity.

Figure 1. Known substances with PAD 4 activity.

Figure 2. Protein Citrullination reaction.

Figure 2. Protein Citrullination reaction.

This study employed two paradigms toward drug designing: structure-based and ligand-based drug design. The study started by using molecular docking to identify the top contenders, which were then further validated by induced-fit docking, molecular mechanics with generalized Born and surface area solvation, and DFT for checking the binding efficiency. Structural compatibility was appraised through shape-based screening. This study has been promising for improving efficacy in treatments directed at PAD4.

2. Materials and Methods

2.1. Computational Tools:

Computational tools: Maestro 13.6 (Schrödinger 2023-2 version) was used to carry out computational simulations [31,32].

2.2. Generation of Databases and Ligand Library Preparation

FDA-approved drugs library was retrieved from the zinc 15 online server (https://zinc.docking.org/), a public web-accessible database containing over 750 million purchasable compounds [33]. 1650 compounds were downloaded and saved in the 2D Structure-Data File (SDF) format and then imported into Maestro. The chirality and ionization states were optimized at a physiological pH of 7.4 ± 2.0 and were determined with the aid of Epik. During this stage, several treatments were applied to the structures. Finally, the OPLS3 force field was selected to optimize the geometries [34].

2.3. Crystal Structure Retrieval and Preparation

The structure of Peptidyl arginine deiminase IV bound to GSK199 was determined using X-ray crystallography and retrieved from the RCSB Protein Data Bank (PDB ID: 4X8G) at a resolution of 3.29 Å as displayed in Figure 3 [25]. The crystal structure was prepared using the Schrödinger Protein Preparation Wizard tool. Adjusting ionization at pH 7.4, adding missing hydrogens, and removing extra water molecules. The protein structure was minimized and optimized by OPLS3 force field to enhance protein energies and avoid steric hindrances, with a default root mean square deviation (RMSD) value of 0.30 Å for non-hydrogen atoms.

2.4. Binding Pocket Determination and Validation of Molecular Docking

In Schrödinger’s Maestro, the binding pocket was identified using the workspace co-crystallized with PAD 4 and GSK 199. The Receptor grid generation tool was then utilized to create a docking grid by selecting the ligand-binding pocket from this crystal structure. A docking grid was established using Glide software, centered on the ligand-binding pocket from the PAD4 co-crystal structure. For validation, the co-crystallized ligand was re-docked again in the same pocket. Superimposition in Maestro and RMSD calculations were used to confirm docking poses and interactions. Receptor grids were generated with a van der Waals radius scaling factor of 1.00 and a partial charge cutoff of 0.25. The docking procedure was subsequently repeated and verified using High Throughput Virtual Screening (HTVS), Standard-Precision (SP) and Extra-Precision (XP) screening settings [35,36,37].

2.5. Standard Molecular Docking (Rigid)

The ligand was docked using the Glide tool without constraints, employing a vdw radius scaling factor of 0.80 and a partial charge cut-off of 0.15. The ligands' flexibility was considered while the protein was considered a rigid structure, with all other parameters set to their default values. GlideScore was utilized to predict ligands' binding affinity. The Pose Rank was utilized to identify the optimal docking pose for each ligand. Binding scores and conformation poses were used in detailed analysis of the resultant structures [32].

2.6. Induced-Fit Docking (IFD) (Flexible)

The IFD technique, developed by Schrödinger, is used to model ligands' binding to different conformational changes. Each ligand undergoes initial docking using a softened potential (van der Waals radii scaling) and flexible conformational sampling. Side-chain prediction is conducted within specified distance of each ligand. Favorable binding poses of structures are predicted based on the IFD score [38].

2.7. Molecular Mechanics-Based Re-Scoring

Molecular mechanics generalized Born surface area (MM/GBSA) docking was used to improve the accuracy of affinity predictions of complexes [39,40].

ΔG binding free energy = ΔG binding, vacuum + ΔG solvation, complex − (ΔG solvation, ligand + ΔG solvation, receptor).

MM/GBSA allow both the ligand and receptor flexibility increasing the accuracy and to relate for normal physiology [41]. An intensive MM/GBSA simulation was used to rank the binding affinities of the five identified hits against the PAD 4 active site. Flexibility was incorporated by adjusting the distance between the hits or GSK 199 and PAD IV to 5 Å. The simulation employed the VGSB solvation model alongside the OPLS3 force field [42].

2.8. Shape-Based Screen

GSK 199 served as the reference structure on Schrödinger's Shape Screening tool. Five compounds underwent screening utilizing the pharmacophore volume scoring technique, which evaluates each compound as an assembly of pharmacophore features, including aromatic groups, hydrogen bond acceptors (HBA), hydrogen bond donors (HBD), hydrophobic regions, as well as positively and negatively charged groups. Shape similarity score was taken from the highest number of matching features [43].

2.9. ADMET Properties and Drug-Likeness Predictions

pkCSM website (http://biosig.unimelb.edu.au/pkcsm/prediction) [44] was used to predict descriptors for ADMET and drug-likeness properties of the final potential inhibitors. Eight molecular descriptors were generated to characterize the ADMET properties of the potential hits. Lipinski's rule of five, the molecular weight, octanol-water partition coefficient (log P), HBDs, and HBAs were applied to predict the physicochemical properties [45].

2.10. Quantum Chemical Calculations

The structural geometries of the four drugs selected with GSK199 for this study were optimized using quantum mechanics density functional theory (DFT) in the ground state with the hybrid functional method B3LYP [46,47], and the LanL2DZ basis set was applied [48]. The Gaussian 09 package was performed to optimize all the structures in the gas phase [49]. The output files of the optimized three-dimensional drugs were visualized using the GaussView 6.0 program [50]. The influential global reactivity descriptors of the chemical structures are determined based on the frontier molecular orbital theory (FMO) consisting of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), as given in equations 1-7. [51,52,53,54,55,56,57]. In addition, the molecular electrostatic potential (MEP) and Mulliken population analysis are scrutinized to investigate the charge distributions and reactivity behaviour of local atoms in selected drugs at B3LYP/ LanL2DZ level of theory. [58,59]

$$\mathbf{1}\mathbf{.}\mathbf{Energy}\mathbf{gab}\mathbf{\left(}\mathbf{Egab}\mathbf{\right)}\mathbf{=}{{\mathit{E}}_{\mathit{L}\mathit{U}\mathit{M}\mathit{O}}\mathbf{-}\mathit{E}}_{\mathit{H}\mathit{O}\mathit{M}\mathit{O}}$$

$$\mathbf{2}\mathbf{.}\mathbf{Ionization}\mathbf{potential}\mathbf{\left(}\mathbf{IP}\mathbf{\right)}\mathbf{=}{\mathbf{-}\mathit{E}}_{\mathit{H}\mathit{O}\mathit{M}\mathit{O}}$$

$$\mathbf{3}\mathbf{.}\mathbf{Electron}\mathbf{affinity}\mathbf{\left(}\mathbf{EA}\mathbf{\right)}\mathbf{=}{\mathbf{-}\mathit{E}}_{\mathit{L}\mathit{U}\mathit{M}\mathit{O}}$$

$$\mathbf{4}\mathbf{.}\mathbf{Global}\mathbf{hardness}\mathbf{(}\mathsf{\eta }\mathbf{)}\mathbf{=}{\displaystyle \frac{\mathit{I}\mathit{P}\mathbf{-}\mathit{E}\mathit{A}}{2}}$$

$$\mathbf{5}\mathbf{.}\mathbf{Global}\mathbf{softness}\mathbf{(}\mathsf{\sigma }\mathbf{)}\mathbf{=}{\displaystyle \frac{\mathbf{1}}{\mathsf{\eta }}}$$

$$\mathbf{6}\mathbf{.}\mathbf{Electronegativity}\mathbf{(}\mathit{c}\mathbf{)}\mathbf{=}{\displaystyle \frac{\mathit{I}\mathit{P}\mathbf{+}\mathit{E}\mathit{A}}{\mathbf{2}}}$$

$$\mathbf{7}\mathbf{.}\mathbf{Electrophilicity}\mathbf{index}\mathbf{(}\mathcal{ώ}\mathbf{)}\mathbf{=}{\displaystyle \frac{{\mathsf{\chi }}^{\mathbf{2}}}{\mathbf{2}\mathsf{\eta }}}$$

3. Results and Discussion

3.1. Docking Studies

The co-crystalline ligand (GSK199) was first redocked into its target PAD4 using the same procedure and protocol applied for the FDA ligands to validate docking. Subsequently, rigid-body superposition was performed using Maestro's structure superposition tool to align the predicted lowest energy conformation of the target with its corresponding co-crystalline ligand. The classical RMSD from the co-crystalline pose was calculated for the predicted binding poses, with an RMSD < 2 Å considered an effective threshold for validating correctly posed molecules [34,35]. The results showed good binding mode superimposition, with an RMSD of 0.527 for GSK199, reflecting the accuracy of Glide's pose prediction Figure 4. XP docking study identified four hits, including Pemetrexed, Leucovorin, Chlordiazepoxide, and Ioversol Figure 5. The binding affinities of the four hits were assessed against PAD4. Initially, the inhibition profiles of these four hits were examined by docking them into the binding pockets of the target, investigating their binding patterns, target interactions, and binding affinities compared to the reference GSK199.

3.2. Computational Analysis of the Four Hits Binding to PAD4

The binding affinities and interactions of the four hit compounds with Peptidyl Arginine Deiminase IV (PAD4) were evaluated using both rigid and induced fit docking (IFD) methodologies (Table 2) [60]. The compounds investigated were Ioversol, Pemetrexed, Chlordiazepoxide, and Leucovorin, with GSK199 serving as a control. Ioversol emerged as the most promising hit, exhibiting superior docking scores and more favorable amino acid interactions than the other tested compounds. The IFD score for Ioversol was -11.617, indicating a robust binding affinity. The compound formed several critical hydrogen bonds with residues within the binding pocket of PAD4, including HID471, ASP473, GLU474, GLU580, ALA581, GLY641, and a conserved water molecule. These interactions are crucial for the stability and specificity of Ioversol within the PAD4 binding pocket. Pemetrexed also demonstrated a high IFD score of -10.599 and established interactions with key residues such as LYS521, LYS572, HIE471, ASP473, ASN585, ASN588, and ALA581. Chlordiazepoxide and Leucovorin displayed IFD scores of -9.988 and -10.521, respectively. Chlordiazepoxide interacted with residues LYS521, PHE633, PHE634, and ALA581, while Leucovorin formed bonds with LYS521, LYS572, LYS521, ALA581, GLU642, PHE633, and HID637. Favorable binding-based 2D and 3D docking positions interact with key residues within the binding pocket, as shown in Figure 6.

3.3. Binding Free Energies Analysis

The MM-GBSA (Molecular Mechanics Generalized Born Surface Area) binding free energy analysis for the four compounds Ioversol, Leucovorin, Chlordiazepoxide, and Pemetrexed, along with the control GSK199, revealed significant differences in their binding affinities to PAD4 (Table 3) [61]. Ioversol has the most favorable net binding free energy of -53.53 kcal/mol, establishing better binding stability and affinity over all other compounds. Leucovorin exhibited a net binding free energy of -43.71 kcal/mol, which was notably less favorable than Ioversol. Chlordiazepoxide and Pemetrexed showed even lower binding affinities, with net binding free energies of -30.96 kcal/mol and -28.46 kcal/mol, respectively.

3.4. Shape Similarity Prediction

Table 4 displays the shape similarity scores of four compounds (Ioversol, Pemetrexed, Chlordiazepoxide, and Leucovorin) that were assessed against GSK199 as a reference [62]. The shape similarity analysis is crucial for understanding how well these compounds mimic the three-dimensional shape of the known inhibitor GSK199, which can influence their ability to bind to the target enzyme, Peptidyl Arginine Deiminase IV (PAD4). Among the compounds evaluated, Pemetrexed demonstrated the highest shape similarity score of 0.300, indicating its structural resemblance to GSK199. Chlordiazepoxide and Leucovorin exhibited shape similarity scores of 0.280 and 0.258, respectively.

3.5. ADMET and Drug-Likeness

The ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) profiling [63] and drug-likeness evaluations show the four compounds under study Ioversol, Pemetrexed, Chlordiazepoxide, and Leucovorin deduce noticeable aspects related to their pharmacokinetic properties and compliance with Lipinski's rule of five (Table 5 and Table 6) [64]. Ioversol exhibited a water solubility of -2.52 and a low Caco2 permeability of 0.17, indicating moderate oral absorption with a human intestinal absorption rate of 22.97%. It showed no significant BBB permeability and CNS permeability, suggesting limited distribution to the brain. Ioversol was not a substrate for CYP3A4 and demonstrated no AMES toxicity, with a maximum tolerated dose of 0.823 log mg/kg/day, and no hERG I inhibition. Pemetrexed showed similar water solubility of -2.87 and Caco2 permeability of -1.16, leading to poor oral absorption with only 15.39% intestinal absorption. It also showed no significant BBB and CNS permeability and was not a CYP3A4 substrate. It was non-toxic in the AMES test, with a maximum tolerated dose of -0.292 log mg/kg/day, and no hERG I inhibition. Chlordiazepoxide which has the least solubility in water at -3.64 had a fair Caco2 permeability of 1.22 thereby resulting in very high human intestinal absorption of 96.32%.

Ioversol did not fully comply with Lipinski's rule of five due to its molecular weight and the number of hydrogen bond donors. Pemetrexed, with a molecular weight of 427.41 and six hydrogen bond donors, had one violation. Leucovorin, with a molecular weight of 473.44 and seven hydrogen bond donors, also had one violation. Chlordiazepoxide fully complied with Lipinski's rule with a molecular weight of 299.76 and acceptable ranges for hydrogen bond donors and acceptors. GSK199, the control, complied fully with the rule, having a molecular weight of 468.98, a log P of 3.60, seven hydrogen bond acceptors, and one hydrogen bond donor (Table 6).

3.6. DFT Optimization Structures

The DFT method is a crucial quantum technique that provides valuable structure optimization details of electronic features relating to the minimum energy of molecules in three-dimensional. This aspect of DFT empowers us with a deeper understanding of the conduct of organic compounds in biological systems, making it an essential tool in many previous studies [65,66,67,68]. Figure 7 shows the optimized geometries of the fourhits and GK199. All compounds were obtained with no imaginary frequency to ensure the drugs were at minimal energy. In addition, the bond distances of all optimized drugs were determined in angstrom units using the B3LYP/LanL2DZ level of theory, as shown in Figure 8.

3.7. Frontier Molecular Orbital

FMO is an essential tool obtained from the DFT quantum calculation method. It describes the energy of the highest orbital consisting electrons (EHOMO) and the energy of the lowest empty orbital level (ELUMO). These energy levels determine a molecule's reactivity and provide substantial evidence regarding its stability. As shown in Figure 11, the HOMO orbital of Ioversol is distributed at the whole molecule except the iodo derivative with side chains of terminals hydroxyl groups. In contrast, the LUMO is attributed to the triiodo phenyl derivatives with amide groups. In addition, the HOMO and LUMO of Chlordiazepoxide were lying around the whole drug except the hydrogen of phenyl groups and LUMO of p-chloro phenyl derivative. Furthermore, the HOMO of Pemetrexed is mainly distributed on 2-amino-3,7-dihydro-4H-pyrrolo[2,3-d] pyrimidine-4-one as the donor group, while the LUMO orbital is found on the acceptor group of the carboxylic chain. Moreover, the HOMO orbital of Leucovorin is mainly located from the formyl piperazine derivative reaching central 4-aminobenzamide moiety. Meanwhile, the LUMO was almost like the LUMO of Pemetrexed. Additionally, the HOMO and LUMO orbitals of GSK199 were in the nearly identical region located from the centre of the Benzimidazole scaffold to the pyrrolo[2,3-b] pyridine terminal.

3.8. Global Chemical Descriptors

The quantum chemical reactivity descriptors are a fundamental approach to scrutinizing the stability and reactivity of molecules according to Koopmans' approximation [69]. The energy gap is a significant parameter that gives information about the stability of molecules. Overall, it is clear from Figure 9 that Chlordiazepoxide had the lowest Egab value (2.892 ev), indicating that it was the most kinetically unstable among the potential drugs, which might be due to the seven-membered ring. Meanwhile, the Egab values of four selected drugs were generally close together. The Egab tendency increases in the following order: Chlordiazepoxide < GSK199 < Pemetrexed < Ioversol < Leucovorin. On the other hand, global hardness and softness are vital criteria for estimating a compound's reactivity and correlate to the Egab. A decline in hardness and a rising softness value indicates a lower energy gap and unstable compounds that are more likely to undergo a reaction. Table 7 shows that Chlordiazepoxide had the lowest hardness and the highest softness value. It was more reactive (unstable) compared to the other selected drugs. In contrast, Leucovorin was the most stable drug, containing a high hardness value and lower softness. At the same time, the GSK199, Pemetrexed, and Ioversol showed moderate reactivity compared to Chlordiazepoxide and Leucovorin. In addition, electronegativity is another helpful parameter that demonstrates the ability of atoms to attract electrons to itself in the chemical bonds in the molecule; the order of the calculated results is as follows, arranged in the growing trends of electronegativity values: Leucovorin < Pemetrexed < GSK199 < Chlordiazepoxide < Ioversol. Moreover, electrophilicity investigates compounds' capability to accept electrons. The electrophilicity index order was diminishing in the rate from Leucovorin> Ioversol > Pemetrexed > GSK199 > Chlordiazepoxide. It was apparent that Ioversol, Pemetrexed, and GSK199 tend to accept electrons from the environment less than Leucovorin and more than Chlordiazepoxide.

3.9. Molecular Electrostatic Potential and Mulliken population Analysis

The molecular electrostatic potential (MEP) is an essential tool for understanding the electron density of molecules around each atom. It is crucial to characterize which atom in a molecule tends to act as an electron donor or acceptor within biological systems. The MEP, as depicted in Figure 10, assists in identifying electron-rich and electron-deficient atoms. The red-to-orange hues represent the negative region, indicating potential electrophilic attack sites due to electron-rich atoms. Conversely, the positive area, indicated by blue to sky-blue cloud, is related to nucleophilic attack regions, representing electron-deficient atoms. Neutral atoms usually appear in white to green shades. Figure 10 provides significant insights into the MEP mapping of selected drugs, the oxygen and nitrogen atoms emerging as the dominant regions of all drugs for electrophilic attack (electron-rich) atoms. At the same time, hydrogens attached directly to highly electronegative atoms are represented as a nucleophilic attack (electron-deficient) region. On the other hand, Mulliken charge is another significant parameter that plays a crucial role in investigating the interaction of molecules with biological systems. It is beneficial in estimating the partial charge distribution at the atom level within a compound. A negative charge indicates an electron-dense atom with a potential region as a strong hydrogen bond acceptor (HBA). Conversely, the atom with positive values efficiently behaves as a hydrogen bond donor (HBD), as shown in Figure 11. The Ioversol, the oxygen of the hydroxyl group (O34), had the highest negative number (-0.537), indicating it had a significant role as an electron-rich atom. Meanwhile, the hydrogen (H35) attached to the same oxygen had the most positive value (0.381). The values of Mulliken atomic charges for the other selected drugs can be found in Figure 11.

Figure 10. 3D of HOMO and LUMO distributions drugs, and the energy level of HOMO, LUMO, and energy gab (ev) values.

Figure 10. 3D of HOMO and LUMO distributions drugs, and the energy level of HOMO, LUMO, and energy gab (ev) values.

Figure 11. Mulliken charges of four optimized geometries drugs with GSK 199.

Figure 11. Mulliken charges of four optimized geometries drugs with GSK 199.

4. Conclusions

The current study presented an in-silico drug repurposing technique of FDA approved drugs for PAD 4 inhibition to identify potential alternative medications for RA. Four hits were identified including Pemetrexed, Leucovorin, Chlordiazepoxide, and Ioversol that showed the best XP scores together with favorable binding interactions. The IFD further supported the results and indicated a robust binding affinity with IFD scores of -11.617, -10.599, -10.521, and -9.988 for Ioversol, Pemetrexed, Leucovorin, and Chlordiazepoxide respectively. The binding free energy calculation proved the favorable binding with scores ranging from -28.46 to -53.53 kcal/mol. The MM-GBSA calculation maintained Ioversol as the first hit while it repositioned Leucovorin, and Chlordiazepoxide as the second and third hits rather than the third and fourth respectively. Furthermore, the binding potential of the four hits and the crystal ligand were analyzed using DFT calculations to ascertain the predicted favorable binding. Of special interest is the identification of Chlordiazepoxide as a potential ligand based on its initial central analgesic activity and tolerability that would position it as a promising candidate for future testing and clinical development.