Molecular interaction of 4-amino-N’-(benzoyloxy)-N-(2,4- dimethylphenyl)-1,2,5-oxadiazole-3-carboximidamide with the methotrexate binding site of human DHFR, and its implication in rheumatoid arthritis

Purpose: To identify an improved lead molecule for the human dihydrofolate reductase (DHFR) inhibition that ‘sits’ in the same binding cavity as methotrexate by high throughput computational screening. Methods: The 3-D structure of the DHFR binding site was examined using ‘CASTp3.0’. Structure based in silico screening of about 5 million drug candidates housed in the MCULE database was performed. The obtained molecule-hits were ranked in accordance with their VINA scores, made to pass through drug-likeness filters, ΔG cut-off criterion, toxicity-checker and finally ‘zero RO5 criterion’. Results: The ‘top molecule’, namely, 4-amino-N'-(benzoyloxy)-N-(2,4-dimethylphenyl)-1,2,5-oxadiazole3-carboximidamide, displayed robust binding with human DHFR through 21 amino acid residues (ΔG = 9.6 kcal/mol) while 10 of these residues were the same as those displayed by ‘methotrexate binding interactions’. It passed through relevant drug screening filters including the ‘Toxicity Checker’. Conclusion: This research work describes the molecular interaction of human DHFR with an improved lead molecule named, 4-amino-N’-(benzoyloxy)-N-(2,4-dimethylphenyl)-1,2,5-oxadiazole-3carboximidamide, with a ΔG of -9.6 kcal/mol, thus satisfying adequate ADME features for further in vitro and in vivo validation in the context of rheumatoid arthritis.


INTRODUCTION
Rheumatoid arthritis (RA) is an autoimmune disorder which predominantly upsets the synovial joints. It leads to disabilities and even early death events. It imposes socioeconomic burden on the ailing individuals as well as on the global community. RA is an area of continued research interest [1] and so is 'Enzoinformatics' [2]. The term 'Enzoinformatics' was coined by Shazi Shakil and accepted by the global scientific fraternity [2]. It is a sub-discipline of Bioinformatics that concerns enzyme-ligand binding interactions in particular [3,4].
The research was driven by the quest of identifying a promising 'seed'-skeleton which could in turn be used for design of potent dihydrofolate reductase (DHFR) inhibitors in near future. It is significant because DHFR-inhibitors have traditionally found application in the treatment of RA. Methotrexate in low dose is the drug of choice in patients with RA, because of its anti-inflammatory action through inhibition of enzymes involved in the folate pathway [12]. Furthermore, application of DHFR-inhibitors for cancer treatment has become an area of increasing research interest [13].
Hence, the objective of our research work was to identify a fresh 'seed-skeleton' for DHFRinhibition that 'sits' in the same binding cavity as methotrexate by high throughput computational screening of about 5 million drug candidates housed in the MCULE database.

CASTp3.0 protein topography probe:
The 3-D structure of the DHFR binding site was examined using 'CASTp3.0' using the method described by Tian et al [14]. It is a free web server. The updated version of this server provides quite intuitive and comprehensive information regarding protein cavities. The protocol in turn utilizes the alpha shape method for identification of relevant protein features, for measurement of volume and area and for computing imprint [14]. The PDB ID 1U72 was duly explored with reference to the precise site of methotrexate binding present on human DHFR using the default probe radius of 1.4Å. As, this one is a complex; the drug was removed by Discovery Studio Visualizer (BIOVIA/Accelerys).
The optimized protein structure equipped for inclusion in the workflow builder was prepared using MAESTRO 9.8 [15].

High throughput computational screening
'Structure-based virtual screening' was performed in MCULE [16] whereby we screened over 5x10 6 putative drug candidates. The objective was to identify a fresh 'seed-skeleton' for DHFR-inhibition that 'sits' in the same binding cavity as methotrexate. A sequential input was fetched to the 'workflow builder' of the MCULE drug discovery platform. In an endeavor to provide initial search flexibility we allowed 1 RO5 violation. A numerical value of 10 was kept for the maximum number of rotatable bonds; minimum value for H-bond donors was entered as 1 while the sampler size was 1000.
An input of 100 was given against 'the maximum number of most-diverse-molecules'. A value of 0.7 was assigned as the threshold similarity cutoff while the rest of the parameters remained as 'default' provided by MCULE. 'Open Babel LF' was employed for the analysis of corresponding descriptors as the screening progressed. The numerical value assigned against "the maximum number of compounds post sphere-exclusion" was 3 x 10 6 as assigned previously [4].

Computational binding experiments
The processed pdb format file (all ligands removed from the complex) of DHFR enzyme was fetched to Vina [17]. The size of the 'Grid' was assigned to be 60 Å × 60 Å ×60 Å 3 for fully covering the methotrexate binding site present on the DHFR protein. The x, y and z grid coordinate values which were essential to mark the 'grid-position' in the three dimensional space were obtained using the information provided with PDB ID 1U72, as done in previous studies for pertinent proteins [4]. Actually, 1U72 is the reality complex describing methotrexate-DHFR interaction. These values for x, y and z as fetched to Vina were 30.620667, 16.876455 and -1.623545, respectively. A default Autodockprotocol was used [17], as done in previous publication [4].

VINA score rankings
The complete set of ligands obtained from the computational screening was ranked as per their VINA scores [17]. It is noteworthy that a ligand which displays a higher (negative) value as its VINA score implies a more efficient binding interaction with the active crevice of the target protein, a principle exploited by several studies [4,5]. Consequently, a list of 50 molecules possessing upper VINA ranks was generated.

SWISS ADME profiling
These 50 ligands were fetched to SWISS ADME to study their comparative drug profiles [18]. This included performing several kinds of tests for the ligands being investigated. Importantly, likelihood of the ligand to behave as a possible future drug molecule was predicted. This likeness to drug was measured by filters, namely Lipinsky (Pfizer), Ghose, Veber (GSK), Egan (Pharmacia) and Muegge (Bayer) filters. The ligands were also tested by 'PAINS' which is an important medicinal chemistry filter. A shortened list of candidate drugs was prepared by rejecting all the ligands that failed against more than 2 filters of drug-likeness.

ΔG cut off
The obtained list of putative drug structures was further shortened by removing all the ligands that displayed a ΔG > -8.6 kcal/mol for their respective complexes with the human DHFR.

Toxicity checker and zero RO5 violation filters
The shortlisted set of candidate drugs was submitted to the 'Toxicity Checker' [16]. This program checks for the presence of possible moieties/sub-structures within the structure of candidate ligands that might be the hall marks of known toxic entities. The toxicity filter passed ligands were further checked for identifying the ligand(s) with zero RO5 violation.

Binding interactions of the 'Top ligand' and 'methotrexate' with DHFR
Discovery Studio 2016 (BIOVIA) was used to detach the methotrexate molecule from its complex with DHFR represented by PDB ID 1U72. Methotrexate was re-docked to the DHFR enzyme by Autodock [17]. The protocol previously published by the corresponding author of this article was also referred [19].
The grid position (x, y and z co-ordinates) was defined by the numerical values as 30.620667, 16.876455 and -1.623545, respectively. A 60 Å 3 grid volume was used for docking. Binding interactions of the 'Top ligand' and 'methotrexate' for their respective complexes with human DHFR enzyme were studied by Discovery Studio Visualizer as well as PyMOL V1.5.0.4 and compared.

High throughput computational screening output
The screening yielded a total of 95 molecules (out of 5 million candidates) in accordance with the provided input in the workflow builder of MCULE ( Figure 1). Fifty molecules that displayed the upper VINA scores were chosen for further filtration. Next, these 50 ligand molecules were tested by Drug-Likeness filters and also the PAINS filter by the SWISS ADME program [18]. Table 1 shows SWISS ADME profiles of the top two inhibitors obtained by screening of 5 million candidate molecules against DHFR enzyme (Table 1).   Table 1]. MCULE-5318527349-0-101, the carboximidamide displayed an XLogP3 value as 3.32. Further increasing the stringency, the two aforementioned ligands were examined for having zero RO5 violations. MCULE-7895890639-0-42 possessed a molecular weight of 510.40 g/mol and hence displayed 1 RO5 violation. Also, the binding energy (negative) was found to be highest for the complex of MCULE-5318527349-0-101 with DHFR (G = -9.6 kcal/mol). Hence, MCULE-5318527349-0-101 i.e. 4-Amino-N'-(benzoyloxy)-N-(2,4-dimethylphe-nyl)-1,2,5-oxadiazole-3-carboximidamide was designated as the 'Top molecule' obtained after the screening in the current study.

Binding interactions of the 'Top ligand' and 'methotrexate' with DHFR:
The interacting amino acid residues for the redocked complex were observed to be the same as that of the reality crystal represented by the PDB ID 1U72. The highest (negative) ΔG value was observed for binding of MCULE-5318527349-0-101 with DHFR (ΔG = -9.6 kcal/mol). The runner up ligand, MCULE-7895890639-0-42 displayed free energy of binding exactly equal to that of the DHFR complex involving methotrexate (ΔG = -8.6 kcal/mol). The poses of these docked complexes, were studied by Discovery Studio Visualizer as well as PyMOL. Figure 2 (Table 2).

DISCUSSION
The human DHFR is a known target for rheumatoid arthritis. Therefore, it is of interest to screen human DHFR for improved lead molecules in comparison with methotrexate [20][21][22]. Such ligand(s) could either act as potent inhibitors themselves or like 'seed-structures available for highly advanced simulation studies. Top VINA scores are generally regarded as indicators of more efficient binding interactions [17].
Ligand molecules which display poor pharmacokinetic profiles tend to be rejected in the initial trail of drug-design process. For imparting initial flexibility to the screening process, at the onset we allowed all ligands that passed at least 4 out of the 6 filters mentioned in Table 1. Also, we allowed 1 RO5 violation for the same reason [23]. A higher (negative) free energy of binding displayed by a docked complex is generally regarded as an indicator of efficient binding interactions for the corresponding protein-inhibitor pair. Therefore, the dissociation rate of such an inhibitor with its binding target would be lower and such inhibitors might be expected to possess an increased half-life [24]. Candidate molecules which possess moieties/sub-structures that might be regarded as signatures or hallmarks of known toxic entities are generally rejected in the virtual screening process. Molecule Leads are normally accepted to have an XLogP3 within 3.5.
Accordingly, MCULE-5318527349-0-101, the carboximidamide displayed an XLogP3 value of 3.32. It exhibited commendable GI-absorption, which is considered as a positive indicator concerning putative drugs that might be given from oral route. It was found to be CNS-inactive, a good feature to avoid off-target effects. The interacting amino acid residues for the re-docked complex were observed to be the same as that of the reality crystal represented by the PDB ID 1U72. Hence, the re-docking results of methotrexate with human DHFR reinforced the accuracy of the dockings performed.
'Amidine' is a more common term used to refer to a 'carboximidamide'. In a 2011 study by Willis et al, it was observed that Cl-amidine was able to partially inhibit arthritis in mice [25]. Rimacalib, although discontinued, was used in trials studying the treatment of Rheumatoid Arthritis and it was a 'carboximidamide' [26]. Amidine has been shown to reduce inflammation and joint destruction in arthritic mice by other authors as well [27]. Methotrexate is understood to have limitations due to toxic effects after long-term use. As DHFR-inhibitors find application in many diseases apart from arthritis, e.g. cancer, significance of novel DHFR-inhibitors or improved lead molecules displaying acceptable toxicity profiles than the contemporary drug candidates is evident.

CONCLUSION
This research work describes the molecular interaction of human DHFR with an improved lead molecule named 4-amino-N'-(benzoyloxy)-N-(2,4-dimethylphenyl)-1,2,5-oxadiazole-3carboximidamide with a ΔG of -9.6 kcal/mol satisfying adequate ADME features for further in vitro and in vivo validation in the context of rheumatoid arthritis.

Acknowledgement
This project was funded by National Plan for Science, Technology and Innovation (MAARIFAH) -King Abdulaziz City for Science and Technology -The Kingdom of Saudi Arabia -award no. 11-MED1553-03. The authors also acknowledge with thanks Science and Technology Unit, King Abdulaziz University for technical support.

Conflict of interest
No conflict of interest is associated with this work.