In silico model and design of novel 5-reductase inhibitors for treatment of benign prostatic hyperplasia

Purpose: To carry out in silico design of 5α-reductase inhibitors and study their potential for use in the treatment of benign prostate hyperplasia (BPH). Methods: In silico molecular docking simulation-based virtual screening of NCI Diversity Set-II containing 1880 diverse ligands was performed against human 5α-reductase for identification of potential lead molecules. The pharmacological properties and toxicity of the lead compounds were determined using software, Marvin Sketch and OSIRIS online programs, respectively. Results: Three compounds: ZINC13099050, ZINC01569237 and ZINC17995347_2 showed potent inhibition of 5α-reductase enzyme protein and good pharmacokinetic properties without any serious toxic effects. Conclusion: The selected lead molecules are promising inhibitors of 5α-reductase. They are recommended for further structure-based development of drugs for the treatment of BPH.


INTRODUCTION
Benign prostate hyperplasia (BPH) is characterised by progressive enlargement of the prostate gland and obstruction of urinary flow by compression of the urethra in elderly men. It results in increased irritability of the detrusor muscle due to chronic obstruction of the urinary tract [1]. The usual symptoms observed in men suffering from BPH are weak urinary stream, urinary urgency, incomplete emptying of the urinary bladder, as well as difficulty in initiating urination, nocturia or dribbling [2]. It was observed in a prevalence survey that more than 42 % of men above 50 years of age have experienced symptoms of BPH, and only one third of this population actually seek treatment for BPH [3,4]. Treatment is not required for patients with mild symptoms of the disease, while the patients with moderate-to-severe symptoms may initially require medical treatment or surgical therapy later [1,2]. Blockage of alpha receptor results in relaxation of the urethral sphincter which helps in emptying the urinary bladder. Vascular smooth muscle relaxation and lowering of blood pressure are enhanced by non-selective alpha blockers, while no effect is observed on resting blood pressure by selective alpha blockers [3,4].
Catheter insertion is required for emptying the urinary bladder in acute urinary retention caused by BPH. Fortunately, only about 1 % men with BPH suffer from the problem of urinary retention annually [1,2]. The use of several categories of medication such as decongestants, general anesthesia, anti-cholinergics, opioids, benzodiazepines, and calcium channel blockers triggers urinary retention in BPH patients [3,4]. Clinical benefits are usually noticeable 2 to 4 weeks after commencement of BPH treatment. The alpha blockers used for the treatment of BPH only help in controlling the symptoms of the disease: they do not affect the pathological basis of the disease since they neither reduce the risk of urinary retention nor the need for surgical intervention [2,3]. A common side effect observed in the use of alpha blockers is orthostatic hypotension. The alpha-blocker drugs are highly contraindicated with drugs used for erectile dysfunction due to increased risk of orthostatic hypotension [1][2][3][4].
Dihydrotestosterone and other androgenic hormones stimulate prostate growth. The conversion of testosterone to dihydrotestosterone is catalysed by 5-alpha reductase, and its inhibition results in reduction in the size of the prostate. The use of 5-alpha reductase inhibitor drugs has resulted in improvement in the control of symptoms of benign prostate hyperplasia. A minimum of six months therapy is required for complete treatment of BPH. Unlike alpha blockers, 5-alpha reductase inhibitors lead to low risk of urinary retention and reduced need for surgical intervention. Thus, there is need for development of newer drugs with selective modes of action and fewer side effects, for the treatment of BPH.

Sequence selection
The primary sequence of human 5α-reductase was retrieved from the NCBI Genbank (Accession Number: AAP36172.1). The 5αreductase enzyme protein sequence was retrieved from Genbank in FASTA format and used for generating its two and threedimensional models [5].

Secondary structure prediction
The secondary structure of 5α-reductase enzyme protein was predicted with FASTA sequences using GOR IV online server. The GOR IV server predicts protein secondary structure from the given protein sequence with respect to annotation of the protein function and its structure [6].

Model building and evaluation
The three-dimensional structure of the enzyme protein was modelled using SWISS MODEL online homology model program. The threedimensional structure model of the protein generated with the SWISSMODEL server depends on the quality of the BLAST sequence alignment and the template structure used. The USFC-Chimera and Swiss PDB Viewer (SPDBV) softwares were used for graphical representation of the protein model, while its evaluation and validation were carried out on the RAMPAGE server through generation of Ramachandran plot [7][8][9].

Prediction of ligand binding site in protein
The ligand binding site of human 5α-reductase enzyme was predicted using SiteHound-web server. The potential ligand binding site was identified with SITEHOUND algorithm through the recognition of certain regions characterized by favourable non-bonded interactions with a chemical probe [10].

Preparation of receptor protein
The tertiary model of 5α-reductase enzyme was prepared for molecular docking simulation by addition of polar hydrogens and Kolhmann's charge, and assigning of AutoDock-4 (AD4) atom types. The format of the tertiary protein model was converted from *.pdb to *.pdbqt using AutoDock [11,18].

Selection of chemical libraries
A total of 1880 diverse molecules present in the NCI "Diversity Set II molecular library were virtually screened to identify possible lead compounds. All the ligands used in virtual screening followed the Lipinsky's rule of five, i.e., ClogP < 5; H bond donor < 5; H bond acceptor < 10 and molecular weight < 500, for good pharmacokinetic properties [12]. Table I shows the grid box points of the receptor protein which were utilized for all docking runs. The grid box was placed in the centre of the ligand-binding site to ensure that all the extended conformations of ligand fitted into it.

Preparation of grid maps
The map files for different atom types in ligands and receptor viz: A, HD, F, I, Br, N, P, S, Cl, SA, NA, C and OA. were prepared by running Autogrid utility of the AutoDock suite. These map files were used in AutoDock for carrying out molecular docking simulations [13,18].

Docking parameters
Lamarckian genetic algorithm (LGA) is the primary conformational search approach in AutoDock [13]. A population was created for various possible conformations of the ligand, and mutations and exchange of conformational parameters of the ligand molecule took place to compete in a manner similar to that of biological evolution in successive generations. This was followed by the final selection of individual ligand with minimum binding energy. The individual conformation searched for local minima and local conformational space. Then, these information were passed to later generations. Semi-empirical free energy force field was used to predict the free energies associated with the binding of small molecules to macromolecular targets. The force field was based on a comprehensive thermodynamic model that allowed intramolecular energy incorporations into the predicted binding free energy by evaluating energies for the bound and unbound states. A docking parameter file was prepared for each ligand using the following parameters: GA runs = 150; maximum number of evaluations (short) = 250000; maximum number of generations = 27000; GA runs 10, and degree of gene mutation = 0.02 [18].

Virtual screening process
The files necessary for performing virtual screening process were prepared using Racoon software. Raccoon is a graphical user interface for AutoDock-based virtual screening of ligand libraries against a macromolecule. It is used to split molecular ligand library file with multiple number of ligands so as to convert them into the pdbqt format required by AutoDock, and filter unwanted ligands using some common criteria such as Lipinski's rules, fragment-like "rule of 3" and drug-likeness. The input file is validated at every step by evaluating for the presence of non-standard atom types and ensuring that parameters, input filenames, and grid maps have coherent formats [13].

Results analysis of virtual screening
After the screening, a script summarize_results.py from Scripps Research Institute was utilized to sort the binding energies of the docked ligands and select the best hits. The results obtained from molecular docking simulation were evaluated on the basis of hydrophobic and polar interactions between ligand and the binding residues present in the active ligand binding site of the macromolecule. Usually, the empirical range of the free binding energy should be in the range of -5 to -15 kcal/mole. The binding affinity was calculated from Eq 1.

K i = e [(ΔG/(RT)] ………… (1)
where ΔG is change in Gibbs free energy upon binding, R is universal gas constant, and T is absolute temperature.

Prediction of physicochemical properties
The various physicochemical properties of the virtually screened ligands such as molecular weight, calculated partition coefficient (ClogP), two-dimensional polar surface area (2D-PSA), number of hydrogen bond donors, and acceptor sites were calculated using software Marvin Sketch [17].

ADME profiling and toxicity prediction of lead compounds
The toxic effects of the virtually screened lead molecules were predicted using OSIRIS online program. The presence of major toxic effects such as mutagenicity, tumorigenicity, irritant effects and reproductive effects in the lead molecules were predicted on the basis of functional groups present in their chemical structures. This program was also used to calculate drug-likeness score and Drug Score of the lead molecules on the basis of their physicochemical properties [14].

Primary structure
The 5α-reductase enzyme protein sequence was retrieved from the NCBI GenBank. The primary protein sequence of 5α-reductase in FASTA format is shown in Figure I.

Secondary structure model
The two-dimensional secondary structure model of 5α-reductase enzyme protein predicted by online server GOR IV is shown in Figure 2. The secondary structure predicted by GOR IV server came out in two outputs forms, one of which is an eye-friendly protein sequence, while the other is in the form of rows present in the predicted secondary structure. The probability values for each secondary structure at each amino acid position are also indicated. The secondary structure was predicted on the basis of the highest probable compatibility with respect to predicted helix segment of a minimum of four residues, and a predicted extended segment of at least two residues.

Tertiary protein model
The tertiary structure of 5α-reductase enzyme protein was predicted using homology model with templates of A chain of protein models 4QUV and 4A2N, while the protein energies of J and N chains of 3J2W model were minimized. The three-dimensional model of 5α-reductase enzyme protein predicted with SWISSMODEL Server is shown in Figure 3.

Protein model evaluation
The three-dimensional model of 5α-reductase enzyme protein was evaluated using RamPage server through generation of Ramachandran plot for the predicted model prepared using SWISSMODEL Server. The Ramachandran plot obtained with RamPage online server displayed all the torsion angles of the protein model. Usually, residue names are displayed on pointing at any location in the plot. This helps to identify the presence of any residues with bad torsion angle in the protein model.

Virtual screening
The potential ligands were selected through analysis of the ligand-protein interactions for top ranking pose of each ligand, and interactions of docked compounds were visualized. The molecular docking results of top five molecules selected in AutoDock-based virtual screening of human 5α-reductase enzyme protein with NCI Diversity set-II ligand library having 1880 small ligand molecules, are shown in Table 2. The best five lead molecules selected on the basis of lowest binding energy after performing molecular docking simulation-based virtual screening of NCI Diversity set-II ligand library were used for determination/calculation of important physicochemical properties, with Marvin Sketch software and Osiris online tool. These properties were partition coefficient (cLogP), two-dimensional polar surface area (2D PSA), molecular weight, as well as hydrogen bond donor (HBD) and hydrogen bond acceptor (HBA) sites. The physicochemical properties of the five best lead molecules for human 5αreductase enzyme protein are shown in Table 3.

ADME and toxicity profiles of lead compounds
The pharmacokinetic parameters of the five lead compounds were determined with respect to absorption, distribution, metabolism, excretion and toxicity, using online server Osiris Molecular Property Explorer. Three of the five selected lead compounds i.e. ZINC13099050, ZINC01569237 and ZINC17995347_2 had good pharmacokinetic profiles, with very little or no toxic effects. The five selected lead molecules were also evaluated for various parameters such as drug-likeness and Drug Score values. The three molecules, i.e., ZINC13099050, ZINC01569237 and ZINC17995347_2 had acceptable levels of druglikeness and Drug Score values. The ADME and toxicity results of the best lead molecules obtained after virtual screening are shown in Figure 5, Figure 6, Figure 7, Figure 8 and Figure  9.   Three out of the five best lead molecules ZINC13099050, ZINC01569237 and ZINC17995347_2 showed promising and potent inhibition of 5α-reductase as well as good pharmacokinetic properties, and they had no toxic effects. Moreover, the three lead molecules had very good drug-likeness and Drug Score values. The ADME and toxicity results of all the lead molecules obtained by virtual screening are shown in Figures 5 -9. The remaining lead molecules showed poor pharmacokinetic profiles, and exhibited some major toxic effects such as mutagenicity, tumorigenicity and adverse reproductive effects.

DISCUSSION
The protein sequence of human 5α-reductase enzyme obtained from NCBI was used to predict its secondary structure on the basis of local information using two networks i.e. sequence-tostructure network and structure-to-structure network. The secondary structure model of the 5α-reductase enzyme was predicted using FASTA sequence with online GOR-IV program. The results revealed that random coil (46.33 %) dominated the secondary structure elements, while alpha helices (27.80 %) and extended strand (25.87 %) were present at nearly equal levels.
The three-dimensional model obtained through homology modelling was evaluated using the Ramachandran plot, which indicated that 84.9 % of residues had φ and ψ angles in the favoured region, while 11.8 % of residues had φ and ψ angles in the allowed regions. The other properties observed in the protein model such as bond angle, bond length and torsion angles were in the range of values expected for a naturally folded protein.
The physicochemical properties of the five best lead molecules for human 5α-reductase enzyme protein showed that three of them had good enzyme inhibition properties and excellent pharmacokinetic profiles, without any major toxic effects such as mutagenicity, tumorigenicity, irritant effects or reproductive effects.

CONCLUSION
A three-dimensional model of the protein sequence of human 5α-reductase has been successfully built with homology modelling technique using online SWISSMODEL server. The predicted model has been validated using Ramachandran plot from RamPage online server. The ligand binding site present in the modelled protein has been identified using SiteHound-web server. The identified ligand binding site was further utilized for the molecular docking simulation-based in silico virtual screening of 1880 diverse ligands in Diversity Set II released by the National Cancer Institute (NCI), using AutoDock. Three compounds ZINC13099050, ZINC01569237 and ZINC17995347_2 show promising results with potent inhibition of 5α-reductase and good pharmacokinetic properties, without the presence of any serious toxic effects. Thus, these molecules are promising lead compounds for further structure-based discovery of novel drugs for the treatment of BPH.