Journal of Biomolecular Structure and Dynamics ISSN: (Print) (Online) Journal homepage: https://www.tandfonline.com/loi/tbsd20 Computational studies on potential small molecule inhibitors of Leishmania pteridine reductase 1 Aaron Boakye, Edward Ntim Gasu, Jehoshaphat Oppong Mensah & Lawrence Sheringham Borquaye To cite this article: Aaron Boakye, Edward Ntim Gasu, Jehoshaphat Oppong Mensah & Lawrence Sheringham Borquaye (2023): Computational studies on potential small molecule inhibitors of Leishmania pteridine reductase 1, Journal of Biomolecular Structure and Dynamics, DOI: 10.1080/07391102.2023.2166119 To link to this article: https://doi.org/10.1080/07391102.2023.2166119 View supplementary material Published online: 12 Jan 2023. Submit your article to this journal View related articles View Crossmark data Full Terms & Conditions of access and use can be found at https://www.tandfonline.com/action/journalInformation?journalCode=tbsd20 JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS https://doi.org/10.1080/07391102.2023.2166119 Computational studies on potential small molecule inhibitors of Leishmania pteridine reductase 1 Aaron Boakyea Borquayea,b , Edward Ntim Gasua,b , Jehoshaphat Oppong Mensaha and Lawrence Sheringham a Department of Chemistry, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana; bCentral Laboratory, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana Communicated by Ramaswamy H. Sarma ABSTRACT ARTICLE HISTORY Leishmaniasis is a neglected tropical disease of major public health concern. Challenges with current therapeutics have led to the exploration of plant medicine for potential antileishmanial agents. Despite the promising activity of some antileishmanial natural products, their protein targets have not been explored. The relevance of folate metabolism in the Leishmania parasite’s existence presents crucial targets for the development of antileishmanial chemotherapy. Pteridine reductase 1 (PTR1), a crucial enzyme involved in DNA biosynthesis, is a validated target of the Leishmania parasite. Unearthing inhibitors of this enzyme is therefore an active research area. The goal of this work is to unearth small molecule inhibitors of PTR1 using molecular docking and molecular dynamic simulations. Thus, the interactions between selected antileishmanial natural products and PTR1 were examined. The binding affinities obtained from molecular docking ranged from 6.2 to 9.8 kcal/mol. When compared to the natural PTR1 substrate biopterin, compounds such as anonaine, chimanine D, corynantheine, grifolin, licochalcone A, piperogalin and xylopine produced better binding affinities, making interactions catalytic residues - Tyr194, Asp181, Phe113, Arg17 and Ser111. The PTR1- xylopine, -piperogalin, -grifolin, and -licochalcone A complexes exhibited remarkable stability under dynamic conditions during the entire 200 ns simulation period. The overall binding free energy of grifolin, piperogalin, and licochalcone A were observed to be 105.711, 103.567, and 105.646 kJ/mol respectively. The binding of these complexes was observed to be favorable and spontaneous and as such capable of inhibiting Leishmania PTR1. They could therefore be considered as candidates in the development of antileishmanial chemotherapy. Received 8 September 2022 Accepted 1 January 2023 Introduction Leishmaniasis is a major public health disease that causes high rates of morbidity and mortality in tropical and subtropical regions around the world (Et-Touys et al., 2017). It is endemic in over 90 countries, including populations in KEYWORDS Leishmaniasis; neglected tropical disease [NTD]; molecular docking; molecular dynamics; grifolin; piperogalin; licochalcone A Africa, Asia, and Latin America. According to the Global Surveillance of Leishmaniasis Report in 2019, over one billion people are at risk in these areas. The disease is known to be spread by roughly 30 Leishmania species including, Leishmania infantum, Leishmania mexicana, Leishmania major, CONTACT Lawrence Sheringham Borquaye lsborquaye.sci@knust.edu.gh Department of Chemistry, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana Supplemental data for this article can be accessed online at https://doi.org/10.1080/07391102.2023.2166119. ß 2023 Informa UK Limited, trading as Taylor & Francis Group 2 A. BOAKYE ET AL. Leishmania brazilienesis, Leishmania aethiopica, Leishmania donovani, and Leishmania amazonensis - using the Phlebotomus sand-fly as the vector. Leishmaniasis exists in three major clinical forms namely visceral (VL), cutaneous (CL), and mucocutaneous leishmaniasis (MCL). Six hundred thousand to one million people are affected by CL and MCL every year (Ruiz-Postigo et al., 2021). VL, which is caused by L. donovani, results in death if left untreated (Gervazoni et al., 2020). Current treatments for leishmaniasis include amphotericin B, pentamidine, and paromomycin, but their toxicity, painful mode of administration, high cost, and issues of resistance have prevented them from being utilized in some endemic areas (Et-Touys et al., 2017; Crentsil et al., 2020; Tulloch et al., 2010). Several plant-based medicines have shown promising activity against the Leishmania parasite. A number of works have successfully isolated the compounds present in these plants and have shown anti-leishmanial activity from in vitro and in vivo screening campaigns. For instance, piperogalin and grifolin from Perperomia galioides have shown activity in vitro against L. donovani, L. braziliensis and L. amazonensis with an IC50 of 10 mg/mL (Mahiou et al., 1995). Corynantheine, dihydrocorynantheine, and ajmalicine, isolated from Corynanthe pachyceras, have also shown in vitro activity against L. major with IC50 values less than 1.5 mM (Staerk et al., 2000) (Table S1). These compounds could potentially be candidate compounds for in silico screening campaigns targeted at unearthing new anti-leishmanial drugs. Despite the significant improvement in finding potential leads for treating leishmaniasis, potential targets of plantderived compounds remain largely unknown. One strategy utilized in the search for lead compounds involves screening a library of compounds against a cellular target of interest. A chosen target should be very crucial for the survival of the parasite. Protein of interest should either be completely absent in the human host, or if present, should possess structural differences that could ensure selectivity for the parasite cell and not the host cell (Hughes et al., 2011). This screening process could be done in vitro as well as in silico. The in silico target-based screening approach have been extensively used to speed up the drug discovery process and was crucial in the design of the first peptidebased HIV protease inhibitors (Monica et al., 2021). This approach could therefore be useful in the search for the next generation of anti-leishmanial drug leads. In the Leishmania parasite, one protein that could potentially be targeted for an in silico target-based screening is pteridine reductase 1 (PTR1) which is very important in the pteridine pathway for DNA synthesis (Tulloch et al., 2010). The pteridine metabolic pathway of the Leishmania parasite is very important in DNA biosynthesis, and is needed for the production of thymidine nucleotides to maintain parasite growth. This pathway provides the largest part of the onecarbon units available to the parasite (Ouellette et al., 2002). The pteridine metabolic pathway is responsible for the reduction of folate and 7,8-dihydrofolate (DHF) to 5,6,7,8-tetrahydrofolate (THF), and the reduction of conjugated and unconjugated pterins into biopterin [DHB] and tetrahydrobiopterin (THB) (Gourley et al., 2001). The THF and THB produced are further converted to N5, N10-methylene tetrahydrofolate by serine hydroxy methyltransferase enzyme (SHMT) through a transamination reaction. The product then serves as a co-factor for thymidylate synthase to catalyze the conversion of deoxyuridine monophosphate to deoxythymidine monophosphate (Tulloch et al., 2010). Folate reductions are catalyzed by dihydrofolate reductase (DHFR) and the pterines are catalyzed by pteridine reductase 1 (PTR1). This suggests that inhibiting enzymes in this pathway will deprive the enzymes of obtaining nucleotides for DNA biosynthesis. Antifolates such as pyrimethamine, sulfa drugs, and trimethoprim have been successful in treating other protozoal infectious diseases but have failed in infections caused by Leishmania, despite this parasite’s extreme divergence from the host (Cavazzuti et al., 2008). Methotrexate has been known to inhibit dihydrofolate reductase (DHFR) but is less active against pteridine reductase 1 (PTR1). PTR1 has been known to serve as a metabolic bypass for the dihydrofolate reductase enzyme. Thus, PTR1 can reduce both folates and pterines in the absence of DHFR. This makes PTR1 essential for the antifolate action (Nare et al., 1997). Gene knock-out studies have highlighted the critical importance of PTR1 in DNA biosynthesis of trypanosomatid protozoans like Leishmania and Trypanosoma (Tulloch et al., 2010; Hardy et al., 1997). Inhibition of PTR1 reduces the cellular pool of deoxythymidine monophosphate, impairing DNA replication and resulting in cell death. The complete absence of PTR1 in the human host presents a prospective protein target for exploration. The active site of PTR1 is an elongated L-shaped cleft of about 22  15 Å, mainly created by C-terminal sections of several b-strands, two sections of a-helix, and an extended loop region (McLuskey et al., 2004). The NADPH cofactor binds in the active site in an extended conformation, in which the nicotinamide ring creates a floor for catalysis, followed by biopterin binding (Nare et al., 1997). Biopterin binds in a manner such that Phe113 dangles just above this catalytic floor, and beneath this Phe113 overhang, lies the pterin-binding pocket (Figure 1). The mechanism for the reduction of biopterin by PTR1 occurs in two phases involving a series of residues. In the first phase reduction, there is a pi-pi stacking between the pterin moiety of the biopterin substrate, the nicotinamide portion of the co-factor, and Phe 113. This orients the substrate such that NADPH is able to transfer a hydride to the C7 atom of biopterin. After hydride transfer, the N8 position of the cofactor is protonated by Tyr 194, which is part of the protein shuffling system (Asp 181 is also part of this system) (Gourley et al., 2001). The first phase reduction leads to the formation of dihydrobiopterin, which serves as a substrate in the second phase reduction. The second phase of reduction of the biopterin proceeds by hydride transfer from the nicotinamide portion of the co-factor to the C6-N5 portion of the dihydrobiopterin formed. The orientation of the dihydrobiopterin is stabilized by hydrogen bonding with Arg 17 and water molecules at the active site. This reaction leads to the formation of tetrahydrobiopterin for further metabolism. The positioning of the cofactor JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS 3 Table 1. Precision docking analysis of selected compounds against PTR1. Compounds Binding affinity Hydrogen bonding interactions Tyr 194[1.85 Å], NAPH 300 [4.01 Å,4.07 Å], Arg 17 [2.83 Å], Ser 111 [2.96 Å] Ser 111[2.71A, 1.97 A], His 38[2.48 Å], Gly 13[2.59 Å] Ser 100 [2.13 Å] Biopterin 8.0 Ajmalicine 9.8 Berberine 8.9 Chimanine D 8.8 Coronaridine 6.5 Corynantheine 8.2 Gly 225[2.08 Å] Dihydrocorynantheine 6.7 Arg 39 [4.27 Å] Dihydrozaluzanin C Diospyrin 6.2 8.4 Grifolin 8.1 Arg 39[1.93 Å] Arg 17[6.74 Å, NAPH 300[3.31 Å,2.24 Å] Asp 181 [2.86 Å] Harmine 8.5 Tyr 171[4.73 Å,2.31 Å,3.70 Å], Asp 158[2.63 Å] Julocrotine 7.9 NAPH 300[2.64 Å], Arg 17[2.40 Å] Licochalcone A 8.3 NAPH 300 [3.65 Å, 4.86 Å] Anonaine 9.5 Muzanzagenin Piperogalin 7.7 9.1 Unonopsine 9.2 Ursolic acid Xylopine 6.5 9.6 Arg 39 [2.01 Å] 2,4,6-triaminoquinazoline [TAQ] 7.8 Asp 181[1.97 Å], NAPH 300[1.83 Å,3.96 Å, 4.10 Å] Asp 142 [2.50 Å] Phe 113[3.60 Å], Tyr 194[2.21 Å], Asp 181[2.41 Å] Phe 113[2.93 Å] Figure 1. Active site cleft of PTR1. The surface of subunit A is shown in grey, the cofactor and inhibitor are represented in stick mode and colored according to atom type: N, blue; O, red; P, pink; C, green for NADPH and red for TAQ. Selected residues that line the active-site cleft are labeled. The image was prepared in the UCSF chimera. Hydrophobic interactions Phe 113[3.94 Å], NAPH 300 [2.17 Å] Asp 142[3.79], Ser 146[2.80A], Leu 18[5.43A]. Pro 201[4.30 Å, 2.56 Å,4.83 Å, Gly 202[3.56 Å] Leu 188[4.35 Å], Tyr 194[5.15 Å], Asp 181[3.12 Å], Phe 113 [3.82,3.08 Å], NAPH 300[3.994 Å] NAPH 300[3.44 Å], Tyr 114[4.15 Å, 4.57 Å], Val 230[4.08, 4.50 Å], Phe 113[4.46 Å], Asp 181[3.48 Å], Tyr 114[3.70 Å,3.91 Å], Asp 142[5.59 Å], Ala 139[5.14 Å]. NAPH 300[5.37 Å] Phe 113[4.68 Å,4.07 Å] Val 230 [5.18 Å, 4.02] NAPH 300[4 Å], Phe 113 [3.65,4.84 Å], Leu 188[5.45 Å,4.14 Å,4.46 Å] Leu 226[4.14 Å], Met 183[5.20 Å] NAPH 261, Met 160[4.41 Å], Gln 163[3.58 Å], Leu 165[2.76 Å4.00 Å], Phe 102 Phe 113[3.84 Å] Met 183[4.86 Å], Leu 188[2.45 Å], Asp 181[2.71 Å], Tyr 194[3.34 Å], Met 183[5.22 Å], Leu 188[3.50 Å, 4.50 Å], Phe 113[4.32 Å], Phe 113[4.57 Å], Asp 181[2.33 Å,2.23 Å], Leu 226, Leu 188[4.933 Å, 3.88 Å], Met 183 [5.32 Å]. Val 230[4.07 Å], NAPH 300[4 Å] Val 230[2.45 Å,5.09 Å] Leu 229, NAPH 300 Asp 181[2.28 Å], Leu 188[3.88 Å,4.93 Å], Met 183, Leu 226[4.64 Å], Phe 113[4.54 Å,3.96 Å] NAPH 300 Phe 113 [3.63 Å, 3.85 Å] involves a series of hydrogen bonds by Ser 111, Lys 198 Val 228, Ser 227 and Gly 225. Lys 198 helps to position the nicotinamide by hydrogen bonding interactions with the cofactor ribose. The basic side chain of Lys 198 may also reduce the pKa of Tyr 194 and in so doing, assist catalysis (McLuskey et al., 2004). One potential requirement needed for PTR1 inhibition is the need for compounds to necessarily interact with these catalytic residues. In effect, the optimum orientations and interactions necessary for the reduction of the biopterin to DHB and THB would be distorted for inhibition of PTR1. In our search for potential anti-leishmanial drug leads, an in silico target based screening approach has been used to identify potential inhibitors of the critical Leishmania protein, PTR1. To do this, natural products with known antileishmanial activities were screened against PTR1 since some of these compounds may probably interact with PTR1 in their antileishmanial action (Figure 2). Molecular docking was used to screen compounds against PTR1. The binding affinities and non-covalent interactions of compounds with PTR1 were 4 A. BOAKYE ET AL. Figure 2. Chemical structures of compounds with antileishmanial activities. evaluated to identify compounds with favorable binding. Selected hits were then subjected to molecular dynamics to interrogate the stability of the binding event in a dynamic environment. We herein report on the binding affinities, complex stabilities, and molecular interactions of antileishmanial natural products against PTR1. The results after molecular dynamics simulation shows favorable binding for grifolin, piperogalin and licochalcone A compared to the biopterin substrate based on their binding energies. This shows that these compounds may bind to the Leishmania PTR1 in order to produce their antileishmanial action. Grifolin, piperogalin and licochalcone A are potential PTR1 inhibitors that could be developed for leishmaniasis chemotherapy. Methods Protein selection and preparation The crystal structure of Leishmania pteridine reductase 1 (PTR1) was retrieved from the Protein Data Bank [https:// www.rcsb.org] with the PDB ID 1W0C and a resolution of 2.60 Å. In this structure, PTR1 exists as a tetrameric structure with 2,4,6-triaminoquinazoline (TAQ), a methotrexate (MTX) mimetic bound in the active site of each monomer (McLuskey et al., 2004). Since the protein is a homotetramer, Chain A was obtained. Protein was prepared by removing crystallographic water and heteroatoms, and keeping the NADPH co-factor in the active site. Using the modeler suit in UCSF Chimera, structures were modeled and refined to replace all missing residues (Sali & Blundell, 1993). Protonation was assigned to histidine, and polar hydrogen atoms were added to rectify the partial charge calculations on standard residues. AMBER ff14SB was used to minimize the energy of the protein backbone structure, while Antechamber in Chimera was used to compute gasteiger charges (Wang et al., 2004; Wang et al., 2001). pKa values were estimated for the unbound protein to identify the ionizable active site residues present on the protein. The PDBID of the protein was inserted into the PROPKA web server (https://www.ddl.unimi.it/vegaol/propka.htm) and the pKa of the various amino acids were obtained (Li et al., 2005). Ligand preparation Structures of ajmalicine, corynantheine, dihydrocorynantheine, coronaridine, anonaine, licochalcone A, piperogalin, grifolin, xylopine, liriodenine, neolitsine, julocrotine, harmine, ursolic acid, muzazangenin, dehydrozaluzanin C, berberine, chimanine D, diospyrin, unonopsine and 2,6-dihydroxy-4methoxychalcone were sketched and modeled using Spartan ’14 (Wavefunction Inc., Irvine California, USA). Geometry optimization and energy minimization were performed using the Hartree-Fork basis set in Spartan ‘14. Compounds were saved in pdb format and were prepared for docking by adding polar hydrogens and gasteiger charges (Borquaye et al., 2020; Mensah et al., 2022). Molecular docking Molecular docking was performed using AutoDock vina (Trott & Olson, 2009) embedded into USCF Chimera. To JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS validate the docking procedure, a co-crystallized ligand, TAQ, was redocked into the active site of PTR1 with the following grid dimensions of X=-5.4661, Y=-5.54676 and Z ¼ 121.997 and size of 21.498 Å. The NADPH cofactor was included in all docking events. The redocked pose and the co-crystallized TAQ were superimposed to determine the extent of deviation. An RMSD value <2 Å confirmed that the procedure was valid for subsequent docking runs. Precision docking was performed with biopterin as the substrate against PTR1, using binding affinity as reference for further analysis. Global docking was performed by setting the grid to cover the entire protein. This was done to determine the preferential binding site of the molecules. Precision docking following the same grid dimensions of TAQ, was then carried out on molecules that were bound at the active site. This was done to exhaustively search and to evaluate the binding affinity of compounds at the active site. After docking, compounds with binding affinities better than biopterin (˂-8.0 kcal/mol) and made key interactions with catalytic residues were selected for molecular dynamic simulation studies. Molecular dynamics simulation study Molecular dynamics (MD) simulations for 200 ns were carried out on grifolin-, piperogalin-, licochalcone A-, anonaine-, xylopine-, corynantheine-, chimanine D- and 2,4,6-triaminoquinazoline-PTR1 complexes as well as apo protein. MD simulations were performed with GROMACS v.2018.6 [19] compiled on the Lengau cluster (Center for HighPerformance Computing, Cape Town, South Africa). The CGenFF (https://cgenff.umaryland.edu/) webserver was used to generate ligand topologies while the protein topology was obtained using the CHARMM36 all-atom force field (Huang & MacKerell, 2013). Solvation was performed using the TIP3P water model (Mark & Nilsson, 2001) in a dodecahedron boundary box. Addition of sodium and chloride ions was done to achieve neutrality. The energies of the systems were minimized to reduce steric clashes using the steepest descent algorithm. Systems were subjected to NVT and NPT equilibration at a constant number, volume, pressure, and temperature ensembles at temperatures up to 300 K and a pressure up to 1 bar for 100 ps. Long-ranged electrostatic interaction was computed using the Particle Mesh Ewald (PME) approach. Cut-offs for Coulomb and van der Waals interactions were set to 1.2 nm. The production run’s defined time steps were 2 fs, with coordinate trajectories written every 10 ps. In all instances, three-dimensional (3D) periodic boundary conditions (PBC) were applied, followed by a 200 ns production run (Kutzner et al., 2015). 5 with VMD 1.9.345 (Visual Molecular Dynamics) using the produced structural (.gro) and pbc compressed trajectory (.xtc) files for all complexes obtained from MD production. Hydrogen bonding analysis All hydrogen bonds, as well as polar hydrogen bonds with bond angles and bond lengths of 180 and 3 Å, respectively, were evaluated to determine the number of hydrogen bonding interactions and their occupancies during the 200 ns MD simulation period. XMGRACE was used to generate the resulting molecular graphics. Binding free MMGBSA) energy calculations (MM–PBSA and The last 50 ns of the protein-ligand complexes were used in computing the average binding energies and the residue-byresidue contribution energies to the overall binding energy using MMPBSA (Kumari et al., 2014) incorporated in GROMACS. The MM-PBSA approach estimates the binding energy by combining (Et-Touys et al., 2017) the energetic terms corresponding to the change in the potential energy in vacuum. This includes bonded terms such as bond, angle, and torsional energies and non-bonded terms such as van der Waals and electrostatic interactions (Ruiz-Postigo et al., 2021), the desolvation of different species, which includes the sum of two energy terms, i.e. the polar and non-polar solvation energy using the implicit solvation model, and (Gervazoni et al., 2020) the configurational entropy associated with the complex formation in the gas phase. Molecular mechanical potential energy (EMM), polar (Gpol), and apolar (Gapol) solvation parameters were included in the estimated binding energies. In a vacuum, EMM estimates each system’s molecular mechanics contribution, which includes bonding (dihedrals, bond stretching, angle bending) and non-bonding (electrostatics and van der Waal’s) interactions. For Gpol, g_ mmpbsa uses the Poisson-Boltzmann Surface Area (PBSA) model to estimate the electrostatic contributions to the solvation free energy, whereas, for Gapol, the solvent-accessible surface area (SASA) design was adopted. These computations were carried out using python scripts (MmPbSaDecomp.py and MmPSaStat.py). The binding energies of the various complexes were estimated as reported by others (Kumari et al., 2014; Kimuda et al., 2019). DGbinding ¼ Gcomplex Gx ¼< EMM > ðGprotein þ Gligand Þ TS þ < Gsolvation > EMM ¼ Ebonded þ Enonbonded ¼ Ebonded þ ðEelectrostatic þ EvdW Þ (1) (2) (3) Gsolvation ¼ Gpolar þ Gnonpolar Stability and trajectory analysis Following the MD productions, the stability of complexes from trajectory analyses were carried out using the simulation output files. Root mean square deviation (RMSD), root mean square fluctuation (RMSF), and radius of gyration (RoG) values were generated and graphed using XMGRACE for stability analysis. Visualizations of trajectories were performed where, G complex is the total free energy of the protein ligand complex and G protein and G ligand are total free energies of the isolated protein and ligand in solvent, respectively. Gx is the protein or ligand or protein ligand complex. hEMMi is the average molecular mechanics potential energy in a vacuum. TS refers to the entropic contribution to the free energy in a vacuum where T and S denote 6 A. BOAKYE ET AL. the temperature and entropy, respectively. The last term hG solvationi is the free energy of solvation. G polar and G nonpolar are the electrostatic and non-electrostatic contributions to the solvation free energy, respectively. To correct for systems that produced a positive binding free energy from g_mmpbsa calculations for stable proteinligand complex(s), the MMGBSA method implemented through the gmx_MMPBSA tool (Vald
es-Tresanco et al., 2021) was used. A single trajectory protocol (STP) was implemented. All energies are reported in kJ/mol. In silico ADME analysis Pharmacokinetic properties of hit compounds from molecular docking including absorption, distribution, metabolism and excretion were determined using the Swiss ADME web server (http://www.swissadme.ch) and ADMET lab 2.0 webserver (https://admetmesh.scbdd.com). Physicochemical parameters depend on factors including molecular weight, number of hydrogen bond acceptor, hydrogen bond donor, number of rotatable bonds, lipophilicity, solubility and topological surface area. Lipinski’s rule of five was used to estimate the medicinal properties of compounds. Absorption parameters include caco2-permeability and human intestinal absorptions. Excretion properties estimated depend on factors such as clearance level and half-life. The SMILES of the compounds were generated from PubChem (https://pubchem.ncbi.nlm.nih.gov) and the pharmacokinetic properties were determined by submitting the smiles to the ADME web server. Results Molecular docking The interactions of the crystallographic pose of TAQ included hydrogen bonding with Tyr 194 (1.58 Å) and Ser 111 (2.37 Å), as well as pi-pi stacking interaction with Phe 113 (3.84 Å, 4.07 Å) (Figure S1a). Docking TAQ to PTR1 gave a binding affinity of 7.8 kcal/mol, following a binding mode which interacted with Asp 181 (1.97 Å) and NAPH 300 (1.83 Å,3.96 Å,4.10 Å) through hydrogen bonding as well as Table 2a. Breakdown of binding energy components obtained from MM-PBSA computations. Protein – Ligand complex Grifolin Piperogalin Licochalcone A Xylopine Biopterin Energies [kJ/mol] DEvdW 137.920 143.933 153.000 137.818 8.371 DEelec 14.409 24.278 23.041 13.503 2.564 DE SOLV 68.724 87.057 94.470 240.245 2.088 DG SASA 22.106 22.413 23.961 21.002 1.307 DG bind 105.711 103.567 105.646 94.929 10.154 DEvdW; van der Waals energy contribution, DEelec; electrostatic energy contribution, DESOLV; Polar solvation energy, DG SASA; Solvent Accessible surface area, DG bind; total binding energy. Phe 113 (3.63 Å,3.85 Å) through pi-pi stacking (Figure S1b). The docked pose of TAQ was superimposed on the co-crystallized TAQ which resulted in a rmsd of 0.016 Å, validating the docking procedure for further applications (Figure S1c). Biopterin, the natural substrate of PTR1, was also docked into the PTR1 active site. A binding affinity of 8.0 kcal/mol, and interactions with Tyr 194 [1.85 Å], NAPH 300 (4.01 Å, 4.07 Å), Arg 17 (2.83 Å), Ser 111 (2.96 Å), Phe 113 (3.94 Å) and NAPH 300 (2.17 Å) were observed. The 21 compounds under investigation were then docked against LmPTR1 using the global docking approach. All compounds were observed to preferentially bind at the active site (Figure S2a). The compounds were therefore subjected to precision docking to assess the binding affinity and interactions with active site residues. The binding affinities of the compounds ranged from 6.2 kcal/mol to 9.8 kcal/mol. Compounds with binding affinities greater than that observed for biopterin (<-8.0 kcal/mol), and having interactions with similar residues such as Tyr 194, Asp 181, Ser 111, Phe 113, and Arg 17 that are known for their catalytic importance and are conserved in PTR 1 orthologs, were selected for further investigation (Table 1). From the results of the docking experiments, anonaine, chimanine D, corynantheine, grifolin, licochalcone A, piperogalin, and xylopine were bound in the active site of PTR1 with binding affinities higher than biopterin (i.e., DGbind ¼ 9.5 kcal/mol, 8.8 kcal/mol, 8.1 kcal/mol, 8.2 kcal/mol, 8.3 kcal/mol, 9.1 kcal/mol and 9.6 kcal/mol respectively). These compounds were also involved in interactions with the highlighted catalytic residues - Tyr 194, Asp 181, Ser 111, Phe 113, and Arg 17. The binding affinities obtained for biopterin and TAQ were 8.0 kcal/mol and 7.8 kcal/mol, respectively. Biopterin interacted through hydrogen bonding with Ser 111, Arg 17, and pi-pi stacking with Phe 113. Similar interactions were also observed in all 7 compounds. Pi-pi stacking interactions with Phe 113 were observed in all 7 compounds. Although interactions with Asp 181 was not observed in biopterin binding, it was observed for grifolin and piperogalin (hydrogen bonding) and in anonaine and xylopine (salt-bridge). Licochalcone A and piperogalin interacted with Tyr 194 through hydrogen bonding. Generally, hydrophobic interactions were observed with Leu 188, Met 183, Leu 226 and Val 230 for all compounds (Figure S2b). The pKa values of the active site residues were generated from the PROPKA web server. Tyr 194, Lys 198, Asp 181 and Arg 17 were observed to be ionizable with pKa values of 15.52, 7.05, 4.91 and 11.75 respectively. Phe 113 and Ser 111 were observed to be non-ionizable (Table S2). Molecular dynamics simulations Anonaine, xylopine, chimanine D, corynantheine, grifolin, piperogalin, licochalcone A, biopterin in complex with PTR1 Table 2b. Breakdown of binding energy components of xylopine bound complex obtained from MM-GBSA computations. Xylopine-bound complex Average VDWAALS 137.82096 EEL 27.02864 EGB 33.22096 ESURF 16.9452 GGAS 110.79232 VDWAALS; van der Waals energy contribution, EEL; electrostatic energy contribution, GSOLV; Polar solvation energy, DG Generalize Born component TOTAL; total binding energy. SURF; GSOLV 50.208 TOTAL 161.00032 non-polar solvation component EGB; JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS as well as PTR1 only (apo protein) were subjected to molecular dynamics simulations. Stability analysis, binding free energy calculations, hydrogen bonding occupancy, residual contribution estimation, and trajectory analysis of protein-ligand complexes were estimated to ascertain the dynamic behavior of the ligands when bound to the protein. Stability analysis The rmsd of the stable ligands relative to their initial docking structural conformation were between 0.5 nm and 1.0 nm, whereas the unstable ligands deviated by over 2 nm. The rmsd for piperogalin-, grifolin-, licochalcone A- and xylopinebound complexes were 1 nm, 0.5 nm, 0.5 nm, and 0.4 nm respectively, relative to the protein backbone obtained from the docking analysis. The rmsd of anonaine, chimanine D and corynantheine were higher (> 2 nm) and were observed to be unstable upon visualization (Figure 3a). The apo protein had a rmsd of 0.5 nm. In the case of all the bound complexes, protein backbone had rmsd ranging between 0.4 nm and 0.6 nm (Figure 3b). Relatively large fluctuations were observed between residues 60-80, 110-145 and 225-250 for both apo protein and 7 ligand-bound complexes, and these regions correspond to loop portions of the protein structure (Figure 2c). Due to their key role in catalysis, the extent of fluctuations of Arg 17, Ser 111, Phe 113, Asp 181, and Tyr 194 were examined throughout the 200 ns simulation (Figure S3). From the residue-specific RMSF analysis of the apo protein, lower fluctuations were observed for Tyr 194. Tyr 194 can be found in the helix region of the protein. Large fluctuations were observed with Arg 17, Ser 111, Phe 113, and Asp 181 when compounds were bound to them. These residues are found on the loop region of the substrate-binding site. For the apo protein, fluctuations were higher in Arg 17 and Phe 113, whereas lower fluctuations of Tyr 194 and Ser 111 were observed. Protein-ligand complexes observed considerable fluctuations compared to apo protein for residues of interest. At residues 60-80 and 110-145, fluctuations were higher for all complexes compared to the apo protein. At residues 225250, fluctuations were lower in xylopine-, TAQ-, chimanine D-, and piperogalin- bound complexes and higher in grifolin-, licochalcone A- and corynantheine-bound complexes. This suggests that the ligand-binding affected the motions of the residues. The radius of gyration of PTR1-grifolin, -piperogalin, -xylopine and -licochalcone A complex ranged from 1.85 nm Figure 3. Stability analysis of protein-ligand complexes. (a) rmsd of ligands relative to their starting conformation. (b) rmsd of protein backbone atoms. (c) RMSF of side chain residues in apo ad protein-ligand complexes. (d) Radius of gyration of apo and bound complexes. 8 A. BOAKYE ET AL. 1.90 nm, and that of the apo protein was 1.80 nm (Figure 3d). Binding free energies Binding free energies for piperogalin-, grifolin-, xylopine-, licochalcone A- and biopterin- bound complexes were calculated using the MMPSA method. The last 50 ns of each simulation was used in the energy computations. The binding free energies of the PTR1-grifolin, -piperogalin, -licochalcone A, -xylopine, and -biopterin complexes were 105.711 kJ/mol, 103.567 kJ/mol, 105.646, 94.929 kJ/mol and 10.154 kJ/mol respectively (Table 2a). From the table, non-polar binding contribution (van der Waals and SASA) to the total binding energies were higher for all compounds. The van der Waals energy contributions of stable complexes, as well as biopterin substrate, ranged from 8.371 to 153 kJ/mol, with the effect greatest in the binding of licochalcone A. High solvation potentials, ranging from 68.724 kJ/mol to 240.245 kJ/mol, were also observed for all ligand-bound complexes. The electrostatic contribution to the total binding energy of the various complexes ranged from 23.041 kJ/mol to 13.503 kJ/mol. GBSA energy calculation was computed for xylopine. Total binding energy of 161.00032 kJ/mol was observed. Van der Waals contribution of 137.82096 kJ/mol was observed to be the highest contributor. Electrostatic, Generalized Born, and non-polar solvation energy contributions were 27.02864 kJ/mol, 7.94 kcal/mol, and 16.9452 kJ/mol respectively (Table 2b). Grifolin with a binding energy of 105.711 kJ/mol showed the most favorable binding under dynamic conditions. For the grifolin-bound complex, interactions with Asp 104 (-3.76 kJ/mol), Thr 165 (-6.48 kJ/mol), Gly 203 (-4.21 kJ/mol), Arg 218 (-4.98 kJ/mol) showed favorable energy contributions whereas an unfavorable binding by Leu 204 (3.03 kJ/mol) was observed. Licochalcone A was stabilized by Arg 102 (-8.36 kJ/mol), Gly 171 (-7.81 kJ/mol), and Leu 210 (-5.93 kJ/mol) with unfavorable contribution from Phe 158 (9.70 kJ/mol). The piperogalin-bound complex was stabilized by Arg 161(-3.78 kJ/mol), Thr 165 (-8.82 kJ/mol), Leu 210 (-6.40 kJ/mol) and Val 257 (-2.64 kJ/mol), with very low unfavorable binding contribution observed for Gln 216 (1.04 kJ/mol). For xylopine, unfavorable binding was observed with Leu 35, Thr 50, Ala 77, Ser 175, and Gln 216 with residue contributions ranging from 8.93 kJ/mol to 21.13 kJ/mol, whereas favorable contributions were observed for Arg 39, Thr 161, Leu 119, Met 179, Gln 22 with contributions ranging from 9.76 kJ/mol to 18.21 kJ/mol (Figure 4). These residues stabilized binding of ligands in the active site allowing compounds to compete with biopterin for binding. The Figure 4. Residue contribution to the total binding energy of protein-ligand complexes. a) Grifolin bound complex b) Licochalcone A bound complex c) Piperogalin bound complex d) Xylopine bound complex. JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS variations in electrostatic contribution were also consistent with polar solvation energy. Hydrogen bonding analysis Due to the contributions of electrostatics to the overall binding energies, hydrogen bonds present in the various proteinligand complexes were evaluated (Table S3). Hydrogen bonding occupancy determines the percentage of simulation time when hydrogen bonds are present. Greater hydrogen bonding interactions increase the stability of protein-ligand complex (Zikri et al., 2020). From the MD simulation, a total of 77 hydrogen bonding interactions were observed in the PTR1-grifolin complex. Highest occupancy occurred between the ligand with Tyr 283, Leu 189 and Lys 244. An average of 2 hydrogen bonds were observed during the simulation. For the PTR1-piperogalin complex, a total of 15 hydrogen bonds were observed. Highest percentage occupancy of 45.07% with Gln 186 was observed. Hydrogen bonding interactions between piperogalin and Asp 181 and His 241 were also observed. The average hydrogen bonds during the entire simulation were 3. In the licochalcone A complex, 79 hydrogen bonds were observed. Highest occupancy of 22.50% with Met 183 was observed. Other hydrogen bonding events with Thr 195, and Trp 283 occurred. An average of 2 hydrogen bonds during the entire simulation were present. Hydrogen bond analysis between PTR1 and xylopine revealed a total hydrogen bond count of 17. Occupancies of 218.64%, 121.94%, and 1.92% were observed with Asp 181, Gly 225, and Ser 227 respectively. Xylopine made an average hydrogen bond of 4 with PTR1 in the simulation period (Figure 5). Analysis of ADME properties The ADME properties of the seven best hit compounds from molecular docking were examined with Swiss ADME and ADMET Lab 2.0 server. Results are explained based on set criteria from the ADMET lab server. The oral bioavailability of molecules was evaluated based on the Lipinski’s rule of five 9 (Molecular weight [MW]  500, lipophilicity [log P]  5, number of hydrogen bond donors [HBDs]  5, and number of hydrogen bond acceptors [HBAs]  10). The absorption parameters included human intestinal absorption, solubility class, andcaco-2-permeability (C2P). Excretion parameters included clearance level, and half-life. If two of the set criteria are violated, poor absorption might be observed. Based on the set conditions, all compounds were observed to obey Lipinski’s rule of five, with moderate to high solubility, and high human intestinal absorption. All compounds produced acceptable C2P values, with high human intestinal absorption properties, and will therefore improve absorption. Poor solubilities, moderate solubilities, and high solubilities were observed in the case of piperogalin and grifolin, licochalcone A and corynantheine, xylopine, and anonaine, respectively. Clearance levels were observed to be higher for piperogalin and grifolin, moderate for xylopine and anonaine, and a low clearance level for licochalcone A and corynantheine. Licochalcone A shows a longer half-life probability. A summary of the pharmacokinetic properties of compounds is shown in Table 3. Discussion Hardy and co-workers highlighted the importance of pteridine reductase 1 in the Leishmania parasite, and thus validated it as a potential target in the folate metabolic pathway for drug candidate exploration (Hardy et al., 1997). Leishmania salvages both folates and unconjugated pterins from the host. The salvage and activation of biopterin and folates occurs by two-successive reduction reactions by hydride transfer. These reactions are facilitated by the pteridine reductase 1 (PTR1) and dihydrofolate reductase thymidylate synthase (DHFR-TS) enzymes. Anti-folate drugs like methotrexates inhibits DHFR-TS but is inactive against PTR1 (Gourley et al., 2001). PTR1 in trypanosomatid protozoans resist the action of methotrexate- a known inhibitor of DHFR. Hence, finding alternative biomolecules that can inhibit PTR1 has become very important. Figure 5. a) Hydrogen bonding in grifolin-, licochalcone a- and piperogalin-bound complexes during simulation. Color code represents each ligand-bound complex as indicated in the legends. b) Hydrogen bonding in xylopine-bound complex. 10 A. BOAKYE ET AL. Table 3. ADME properties of best hit compounds from molecular docking study. Physicochemical Compounds MW Grifolin Piperogalin Chimanine D Licochalcone A Xylopine Corynantheine Anonaine 328.49 328.24 185.22 338.15 295.33 366.19 265.31 n HA n HD 2 2 2 4 4 4 3 2 2 0 2 1 1 1 Medicinal n Rot cLog P TPSA Lipinski 8 7 1 6 1 5 0 5.82 5.83 2.39 4.205 2.88 3.08 2.88 40.46 40.46 25.42 66.76 39.72 54.56 30.49 Accepted Accepted Accepted Accepted Accepted Accepted Accepted Absorption range C2P 4.736 4.738 4.588 4.862 4.917 4.43 4.791 HIA High High High High High High High Excretion Clearance T1/2 14.88 17.096 2.529 2.904 10.122 3.997 12.045 0.199 0.221 0.186 0.741 0.129 0.229 0.115 Solubility ESOL Log 6.08 6.19 2.79 4.98 3.77 4.3 3.71 Solubility Class Poorly soluble Poorly soluble Soluble moderately soluble Soluble moderately soluble Soluble MW; Molecular weight, nHA; number of hydrogen bond acceptors, nHD; number of hydrogen bond donor, nRot; number of rotatable bonds, cLog P; lipophilicity, C2P; caco-2-permeability, HIA; Human Intestinal Absorption, T1/2; Half-life. Molecular docking and molecular dynamics simulations have been utilized in finding potential inhibitors of PTR1. A team recently screened non-nucleoside compounds from the ZINC database against LmPTR1. The approach identified Z80393 as potential inhibitor whose inhibition mechanism was investigated by thermal shift assay. Molecular docking and molecular dynamics were used to analyze the interaction profiles of Z80393 (Leite et al., 2017). Kimuda and co-workers also performed a molecular docking of a ligand library of 5742 compounds against TbPTR1 and LmPTR1. Eighteen of the compounds showed promising binding modes and were then subjected to molecular dynamics to characterize their molecular interactions and energetics, followed by in vitro testing. From the results obtained, ZINC00809143 [N’-f[1-[2,4-dimethylphenyl]-3-methyl-5-oxo-1,5-dihydro-4H-pyrazol-4ylidene]methylg-3-nitrobenzohydrazide], ZINC00630525 [N-[4methoxyphenyl]-2-[[4-oxo-5-phenyl-4,5-dihydro-1H-pyrazolo[3,4d]pyrimidin-6-yl]sulfanyl]acetamide], ZINC0058117 (eriodictyol) and ZINC04313814 [2-[2-[cyclopentylidenehydrazono]-4hydroxy-2,5-dihydro-1,3-thiazol-5-yl]-N-phenylacetamide] showed promising activities both in vitro and in silico (Kimuda et al., 2019). These observations highlight the importance of molecular docking and molecular dynamics in unearthing potential drug candidates. In the reduction of biopterin by PTR1, the importance of Arg 17, Ser 111, Phe 113, Asp 181, Tyr 194, Lys 198, and Arg 287 in creating the active site and also facilitating substrate binding has been emphasized (Gourley et al., 2001; Nare et al., 1997; Hardy et al., 1997). Crystallographic analysis of the binding of TAQ shows that TAQ mimics the pterin fragment of methotrexate and inhibit PTR1 with an IC50 comparable to methotrexate (Hardy et al., 1997; McLuskey et al., 2004). TAQ and biopterin interact with catalytic residues such as Tyr 194 and Phe 113. These observations suggest that potential inhibitors must interact with these residues or show good binding affinity than biopterin. When this happens, the optimal orientation of the substrate needed for catalysis to proceed will be distorted. The docking of biopterin to PTR1 proceeded with a binding affinity of 8.0 kcal/mol. Compounds with good binding affinities greater than 8.0 kcal/mol could potentially act as a competitive inhibitor. Based on these, binding affinity and interacting residues are critical factors in identifying potential inhibitors of PTR1. The binding affinities of the molecules evaluated in the molecular docking study ranged from 6.2 kcal/mol to 9.8 kcal/mol. In general, higher binding affinities were observed for the aporphine and indole alkaloids, and those that possessed a phenolic scaffold. Compounds with the aporphine scaffold possessed the highest binding affinities, likely because of the extensive hydrophobic interactions the framework made with the protein. Differences in binding affinities of compounds with the aporphine backbone originated largely from substitutions on the backbone (Table 1). Hardy et. al., reported that the quinazoline moiety can mimic the pterin moiety of biopterin in the PTR1 active site (Hardy et al., 1997). The greater binding affinity of aporphine alkaloids may likely be due to functional group similarities between aporphine alkaloids and the quinazoline structure. It has been reported that binding affinity increases with an increase in hydrophobic interactions (Mohapatra et al., 2015). The compounds with the phenolic backbone also interacted with PTR1 with very good binding affinities. The aromatic group was crucial for the pi-pi stacking interaction with Phe 113. Carbonyl and alcohol functionalities mediated hydrogen bonding interactions with the catalytic residues. Aporphine alkaloids (anonaine and xylopine) phenolics (grifolin, piperogalin and licochalcone A) chimanine D, and corynantheine (indole alkaloid) all exhibited very good binding affinities compared to biopterin. These compounds interacted with Tyr 194, Asp 181, Phe 113, Ser 111, and other hydrophobic residues such as Leu 188, Val 230, Leu 266, and Met 183, and hence were selected for further investigations. The probable protonation states of the active site residues were analyzed using the PROPKA web server. From the results, Tyr 194 had a pKa of 15.52. Therefore, the ionizable proton is easily donated [39]. Interestingly, Tyr 194 donates a hydrogen atom to Asp 181 for biopterin reduction mechanism (Gourley et al., 2001). These observations however correlate with experimental prediction. Lys 198, Asp 181 and Arg 17 had positive pKa values. This suggest that protons present are not easily donated. Although molecular docking provides the lowest energy and binding conformation of the ligand bound to the protein, it does not consider the conformational changes of the protein, solvation, and ionic effects that occur in nature, as well as the entropic contributions to the total binding free energies (Leite et al., 2017; Kyei et al., 2022). Molecular dynamics accounts for these limitations in docking and estimates the stability and corresponding energies of ligands when bound to their targets. Stability analysis, binding free energy calculations, hydrogen bonding occupancy, residual contribution estimation, and trajectory analysis of protein-ligand complexes were estimated to interrogate the dynamic JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS 11 Figure 6. Binding interactions of xylopine (a) and anonaine (b) after temperature and pressure equilibrations. Frame extracted at 0 ns of the simulation time. behavior of the ligands when bound to the protein (Karplus & Kuriyan, 2005). Conformational stability of protein-ligand complexes was estimated using their root mean square deviation (rmsd), root mean square fluctuation [rmsf], and radius of gyration (Rg). The radius of gyration determines the overall compactness of the protein. A low Rg is an indication that the protein does not undergo significant fluctuations, therefore resulting in higher compactness (Justino et al., 2021). The Rg of all protein-ligand complexes did not deviate significantly in comparison to the apo protein (Figure 3d). This suggests that ligand binding did not disturb the overall fold of PTR1. These observations are consistent with experimental data presented by Hardy. PTR1 does not undergo a large conformation change upon ligand binding due to the rigidity of the active site (Hardy et al., 1997). The rmsf measures the extent of side chain deviations of amino acid residues (Mart
ınez, 2015). Similar regions of protein fluctuations were observed for both apo protein and the ligand bound complexes. Fluctuations were observed within residues 60 to 80, 110 to 140, and around 225 to 250 (Figure 3c). These regions correspond to loop regions in the protein. Loop regions undergo conformational changes in most catalytic enzymes, and significantly shields the active site from the solvent. This alters the chemical environment around the substrate (Hospital et al., 2015). The rmsf and Rg data shows that the active site of PTR1 is rigid and ligand binding do not significantly affect fluctuations of active site residues. This correlates with the active site rigidity reported by Hardy et al. In the rmsd plot of the protein backbone of the complex relative to the apo protein, deviations compared to the apo protein were not significant due to the rigidity of the protein. The rsmd of the co-crystallized inhibitor, TAQ was >2 nm, and hence unstable. This was confirmed by visualizing interactions between TAQ and PTR1 after 10 ns of the simulation time. It’s been reported that the multi-targeting nature of the 4-aminopyrimidine derivatives makes them poor affinity binders (Patil et al., 2010). This may be the reason TAQ does not bind to PTR1 during the entire simulation. Analyses of the PTR1-anonaine complex reveals that anonaine was initially bound to PTR1, only to be kicked out after 80 ns of the simulation. Interestingly, the PTR1-xylopine complex was relatively stable throughout the simulation even though xylopine shares structural similarity to anonaine. Comparing both structures after temperature and pressure equilibration, both xylopine and anonaine established interactions with Leu 188, Phe 113 and Tyr 194. Xylopine, however, made extra interactions with Tyr 283, His 241 and Leu 226, and these interactions contributed to the stability of the compound (Figure 6). Compounds with the phenolic scaffold (piperogalin, grifolin, and licochalcone A) were also observed to have a stable rmsd plot as well. The aromatic ring and alcohol functionalities present on the phenolic scaffolds were very important for the binding even under dynamic conditions. A lower rmsd of licochalcone A was observed because the double bond restricted the rotation of the two benzene rings. Binding free energies were estimated for the stable complexes. Using the last 50 ns, the binding energies of grifolin-, licochalcone A-, piperogalin-, xylopine-, and biopterin-bound complexes were computed. The total binding free energy of the individual complexes is as a result of the individual electrostatic and van der Waals interactions between protein and ligand. Solvation contributes unfavorably to ligand binding, and could enhance the dissociation of bound complexes. Grifolin-, piperogalin-, xylopine, licochalcone A, and biopterin-bound complexes showed binding energies of 105.711 kJ/mol, 103.567 kJ/mol, 94.929 kJ/mol, 105.646 kJ/mol and 10.154 kJ/mol. Binding was favorable and spontaneous in grifolin-, piperogalin- and licochalcone A-bound complexes. Grifolin, licochalcone A and piperogalin showed better binding energies than biopterin. Positive binding free energy was observed for xylopine. Due to structural rigidity of xylopine-PTR1 active site complex, low conformational changes were observed. This resulted from very strong interactions with residues, which affected the 12 A. BOAKYE ET AL. conformational (rotational and vibrational) degrees of freedom of protein side chains contributing to the overall low entropy (Bronowska, 2011). Kyei et. al. also reported a positive binding free energy for cryptoheptine and cryptolepinone against Plasmodium falciparum DHFR. They attributed their observation to the inability of MMPBSA to sample varying conformations of ligands to account for their entropic contributions (Kyei et al., 2022). MM-GBSA calculation was employed for xylopine-bound complex. Binding energy of 161.00032 kJ/mol was observed (Table 2b). This suggests that binding event was favorable as opposed by MMPBSA result. The van der Waals (DEvdw), the nonpolar component of the solvation energy (DEsurf) and the polar component (DEGB) of solvation free energy (DGsolv), terms of the molecular mechanics energy showed favorable complex stability and ligand interaction. The most important contributions to the total binding energy are estimated from the van der Waals energy and electrostatic energy. Electrostatic contribution arises from the hydrogen bonding interactions directly from the ligands or due to the presence of water-mediated hydrogen bonds (Kumari et al., 2014). Due to the electrostatic contributions, the extent of hydrogen bonding of the piperogalin-, grifolin-, and licochalcone A-bound complexes were estimated. Interesting trends were observed between electrostatic contributions and hydrogen bonding interactions for the phenolic scaffolds. Piperogalin, with the highest electrostatic contribution, produced the highest average number of hydrogen bonds. Grifolin had the lowest electrostatic contributions and produced the lowest average number of hydrogen bonds. In all complexes, van der Waals energy were significant. This suggests the importance of hydrophobic interactions. The DGSASA (Solvent Accessible Surface Area) analysis was performed to estimate the non-electrostatic term of the polar solvation energy. These include the repulsive and attractive forces between the solute and the solvent, generated by the van der Waals interactions (Kumari et al., 2014). SASA energies were similar in all cases and had a positive effect on overall binding energy. The overall binding energy components explain the importance of hydrophobic core, aromatic rings, and hydrogen bond donor-acceptor groups in ligand binding. Due to the structural differences between piperogalin, grifolin, licochalcone A, and xylopine from biopterin, residual contributions to the total binding energies of ligands were also explored. The ADME properties of anonaine, xylopine, grifolin, piperogalin, licochalcone A, chimanine D, and corynantheine were analyzed to evaluate their pharmacokinetic properties (Table 3). Based on Lipinski’s rule of 5, human intestinal absorption, caco-2-permeability, clearance level, solubility class, and halflife. Computational approach to evaluate ADME, and toxicity parameters of compounds is important in drug discovery processes (Boobis et al., 2002). According to the rule of 5, a molecule would not be orally active if it violates two or more of the four rules (Guan et al., 2019). All compounds did not violate any of Lipinski’s rule and are therefore likely to be orally bioavailable. Absorption measures movement of a drug from an extravascular site of administration into the systemic circulation (Dur
an-Iturbide et al., 2020). Absorption depends on the solubility, physicochemical properties, and lipophilicity of the compound. Insufficient solubility may limit the Human Intestinal Absorption (HIA) through the portal vein system to obtain a therapeutic effect when systemic effects are warranted (Dur
an-Iturbide et al., 2020; Daina et al., 2017). High HIA suggest that compounds can be orally absorbed by the intestine. Caco-2-permeability access the rate at which a compound accesses a fully differentiated monolayer of human epithelial cells (Stockdale et al., 2019). Acceptable HIA and caco-2 properties indicate sufficient solubility, physicochemical property, and lipophilicity of compounds, hence acceptable absorption properties. Excretion of compounds influence both the half-life and bioavailability thus impacting the dose regimen and dose size of a drug (Dur
an-Iturbide et al., 2020). The low clearance level observed for licochalcone A suggests that it can stay for a longer time for absorption to occur. Further optimizations are recommended to improve the solubilities of grifolin and piperogalin. Conclusion This work implemented molecular docking and molecular dynamics simulation studies to screen anti-leishmanial natural products as potential inhibitors of Leishmania pteridine reductase 1 (PTR1). Seven compounds exhibited strong binding affinities against PTR1 and were therefore subjected to molecular dynamics simulations. Four of the seven compounds were stable in the active site under dynamic conditions. The MMPBSA calculations suggested spontaneous binding of three compounds with PTR1, further characterizing them as potential competitive inhibitors. The drug-like properties of the three compounds were highlighted by ADME analysis, and their research as potential antileishmanial drugs is strongly recommended. Acknowledgement The authors are grateful to the Center for High Performance Computing, Cape Town, South Africa who granted generous access to the Lengau cluster for the MD simulations. Disclosure statement All authors declare that there were no financial, professional, or personal competing interests that might have influenced the performance or presentation of this study. Funding The author(s) reported there is no funding associated with the work featured in this article. ORCID Aaron Boakye http://orcid.org/0000-0001-6040-8671 Edward Ntim Gasu http://orcid.org/0000-0002-1752-0948 http://orcid.org/0000-0003-4676-9299 Jehoshaphat Oppong Mensah http://orcid.org/0000-0002-5037-0777 Lawrence Sheringham Borquaye JOURNAL OF BIOMOLECULAR STRUCTURE AND DYNAMICS Author contributions LSB, AB and ENG conceived the study. All experiments were designed by LSB, AB, JOM and ENG. Computations were made by AB, JOM and ENG. Data analysis was by AB, ENG, JOM and LSB. Manuscript was prepared by AB, and edited by all the authors. All authors read and approved the final manuscript. Data availability statement All data generated or analyzed during this study are included in this published article and its supplementary materials. References Boobis, A., Gundert-Remy, U., Kremers, P., Macheras, P., & Pelkonen, O. (2002). In silico prediction of ADME and pharmacokinetics: Report of an expert meeting organised by COST B15. European Journal of Pharmaceutical Sciences: Official Journal of the European Federation for Pharmaceutical Sciences, 17(4-5), 183–193. Borquaye, L. S., Gasu, E. N., Ampomah, G. B., Kyei, L. K., Amarh, M. A., & Mensah, C. N. (2020). Alkaloids from Cryptolepis sanguinolenta as potential inhibitors of SARS-CoV-2 viral Proteins: An in silico study. Biomed Research International. 2020, e5324560. Bronowska, A. K. (2011). Thermodynamics of ligand-protein interactions: implications for molecular design. In Thermodynamics-interaction studies-solids, liquids and gases. IntechOpen. Cavazzuti, A., Paglietti, G., Hunter, W. N., Gamarro, F., Piras, S., Loriga, M., Allecca, S., Corona, P., McLuskey, K., Tulloch, L., Gibellini, F., Ferrari, S., & Costi, M. P. (2008). Discovery of potent pteridine reductase inhibitors to guide antiparasite drug development. Proceedings of the National Academy of Sciences of the United States of America, 105(5), 1448–1453. Crentsil, J. A., Yamthe, L. R. T., Anibea, B. Z., Broni, E., Kwofie, S. K., Tetteh, J. K. A., & Osei-Safo, D. (2020). Leishmanicidal potential of hardwickiic acid isolated from croton sylvaticus. Frontiers in Pharmacology, 11, 753. Daina, A., Michielin, O., & Zoete, V. (2017). SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules. Scientific Reports, 7(1), 42717. https:// doi.org/10.1038/srep42717 Dur
an-Iturbide, N. A., D
ıaz-Eufracio, B. I., & Medina-Franco, J. L. (2020). In silico ADME/Tox profiling of natural products: A focus on BIOFACQUIM. ACS Omega, 5(26), 16076–16084. https://doi.org/10. 1021/acsomega.0c01581 Et-Touys, A., Bouyahya, A., Fellah, H., Mniouil, M., El Boury, H., Dakka, N., Sadak, A., & Bakri, Y. (2017). Antileishmanial activity of medicinal plants from Africa: a review. Asian Pacific Journal of Tropical Disease, 7(12), 826–840. https://doi.org/10.12980/apjtd.7.2017D7-215 Gervazoni, L. F. O., Barcellos, G. B., Ferreira-Paes, T., & Almeida-Amaral, E. E. (2020). Use of natural products in Leishmaniasis chemotherapy: An overview. Frontiers in Chemistry, 8, 579891. https://www.frontiersin. org/articles/10.3389/fchem.2020.579891/full Gourley, D. G., Sch€ uttelkopf, A. W., Leonard, G. A., Luba, J., Hardy, L. W., Beverley, S. M., & Hunter, W. N. (2001). Pteridine reductase mechanism correlates pterin metabolism with drug resistance in trypanosomatid parasites. Nature Structural Biology, 8(6), 521–525. Guan, L., Yang, H., Cai, Y., Sun, L., Di, P., Li, W., Liu, G., & Tang, Y. (2019). ADMET-score – A comprehensive scoring function for evaluation of chemical drug-likeness. MedChemComm, 10(1), 148–157. Hardy, L. W., Matthews, W., Nare, B., & Beverley, S. M. (1997). Biochemical and genetic tests for inhibitors of Leishmania pteridine pathways. Experimental Parasitology, 87(3), 158–170. ~i, J. R., Orozco, M., & Gelp
ı, J. L. (2015). Molecular Hospital, A., Gon dynamics simulations: Advances and applications. Advances and Applications in Bioinformatics and Chemistry, 8, 37. Huang, J., & MacKerell, A. D. Jr (2013). CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data. Journal 13 of Computational Chemistry, 34(25), 2135–2145. https://doi.org/10. 1002/jcc.23354 Hughes, J. P., Rees, S., Kalindjian, S. B., & Philpott, K. L. (2011). Principles of early drug discovery. British Journal of Pharmacology, 162(6), 1239– 1249. https://doi.org/10.1111/j.1476-5381.2010.01127.x Justino, G. C., Nascimento, C. P., & Justino, M. C. (2021). Molecular dynamics simulations and analysis for bioinformatics undergraduate students. Biochemistry and Molecular Biology Education: A Bimonthly Publication of the International Union of Biochemistry and Molecular Biology, 49(4), 570–582. https://doi.org/10.1002/bmb.21512 Karplus, M., & Kuriyan, J. (2005). Molecular dynamics and protein function. Proceedings of the National Academy of Sciences of the United States of America, 102(19), 6679–6685. https://doi.org/10.1073/pnas. 0408930102 € (2019). Kimuda, M. P., Laming, D., Hoppe, H. C., & Tastan Bishop, O. Identification of novel potential inhibitors of pteridine reductase 1 in Trypanosoma brucei via computational structure-based approaches and in vitro inhibition assays. Molecules, 24(1), 142. https://doi.org/10. 3390/molecules24010142 Kumari, R., Kumar, R., & Lynn, A. (2014). g_mmpbsa—A GROMACS tool for high-throughput MM-PBSA calculations. Journal of Chemical Information and Modeling, 54(7), 1951–1962. Kutzner, C., P
all, S., Fechner, M., Esztermann, A., de Groot, B. L., & Grubm€ uller, H. (2015). Best bang for your buck: GPU nodes for GROMACS biomolecular simulations. Journal of Computational Chemistry, 36(26), 1990–2008. https://doi.org/10.1002/jcc.24030 Kyei, L. K., Gasu, E. N., Ampomah, G. B., Mensah, J. O., & Borquaye, L. S. (2022). An in silico study of the interactions of alkaloids from Cryptolepis sanguinolenta with Plasmodium falciparum dihydrofolate reductase and dihydroorotate dehydrogenase. Journal of Chemistry. 2022, e5314179–26. https://doi.org/10.1155/2022/5314179 Leite, F. H. A., Froes, T. Q., da Silva, S. G., de Souza, E. I. M., Vital-Fujii, D. G., Trossini, G. H. G., Pita, S. S. d R., & Castilho, M. S. (2017). An integrated approach towards the discovery of novel non-nucleoside Leishmania major pteridine reductase 1 inhibitors. European Journal of Medicinal Chemistry, 132, 322–332. Li, H., Robertson, A. D., & Jensen, J. H. (2005). Very fast empirical prediction and rationalization of protein pKa values. Proteins: Structure, Function, and Bioinformatics, 61(4), 704–721. https://doi.org/10.1002/ prot.20660 Mahiou, V., Roblot, F., Hocquemiller, R., Cav
e, A., Barrios, A. A., & Fournet, A. (1995). Piperogalin, a new prenylated diphenol from Peperomia galioides. Journal of Natural Products, 58(2), 324–328. Mart
ınez, L. (2015). Automatic identification of mobile and rigid substructures in molecular dynamics simulations and fractional structural fluctuation analysis. PloS One, 10(3), e0119264. McLuskey, K., Gibellini, F., Carvalho, P., Avery, M. A., & Hunter, W. N. (2004). Inhibition of Leishmania major pteridine reductase by 2,4,6-triaminoquinazoline: Structure of the NADPH ternary complex. Acta Crystallographica. Section D, Biological Crystallography, 60(Pt 10), 1780–1785. Mensah, J. O., Ampomah, G. B., Gasu, E. N., Adomako, A. K., Menkah, E. S., & Borquaye, L. S. (2022). Allosteric modulation of the main protease (MPro) of SARS-CoV-2 by casticin—insights from molecular dynamics simulations. Chemistry Africa, 5(5), 1305-1320. https://doi. org/10.1007/s42250-022-00411-7 Mohapatra, S., Prasad, A., Haque, F., Ray, S., De, B., & Ray, S. S. (2015). In silico investigation of black tea components on a-amylase, a-glucosidase and lipase. Journal of Applied Pharmaceutical Science 5(12), 042– 047. Monica, L., Lauria, G., Bono, A., & Martorana, A. (2021). A. off-targetbased design of selective HIV-1 PROTEASE inhibitors. International Journal of Molecular Sciences, 22(11), 6070. https://doi.org/10.3390/ ijms22116070 Nare, B., Hardy, L. W., & Beverley, S. M. (1997). The roles of pteridine reductase 1 and dihydrofolate reductase-thymidylate synthase in pteridine metabolism in the protozoan parasite Leishmania major. The Journal of Biological Chemistry, 272(21), 13883–13891. https://doi.org/ 10.1074/jbc.272.21.13883 14 A. BOAKYE ET AL. Ouellette, M., Drummelsmith, J., El Fadili, A., K€ undig, C., Richard, D., & Roy, G. (2002). Pterin transport and metabolism in Leishmania and related trypanosomatid parasites. International Journal for Parasitology, 32(4), 385–398. https://doi.org/10.1016/s00207519(01)00346-0 Patil, R., Das, S., Stanley, A., Yadav, L., Sudhakar, A., & Varma, A. K. (2010). Optimized hydrophobic interactions and hydrogen bonding at the target-ligand interface leads the pathways of drug-designing. PloS One, 5(8), e12029. https://doi.org/10.1371/journal.pone.0012029 Ruiz-Postigo, J. A., Jain, S., Maia-Elkhoury, A. M. A. N., Valadas, S., Warusavithana, S., Osman, M., et al. (2021). Global Leishmaniasis surveillance: 2019–2020, a baseline for the 2030 roadmap/Surveillance mondiale de la leishmaniose: 2019–2020, une periode de reference pour la feuille de route a l’horizon 2030. Weekly Epidemiological Record, 96(35), 401–420. Sali, A., & Blundell, T. L. (1993). Comparative protein modelling by satisfaction of spatial restraints. Journal of Molecular Biology, 234(3), 779– 815. Staerk, D., Lemmich, E., Christensen, J., Kharazmi, A., Olsen, C. E., & Jaroszewski, J. W. (2000). Leishmanicidal, antiplasmodial and cytotoxic activity of indole alkaloids from Corynanthe pachyceras. Planta Medica, 66(6), 531–536. Stockdale, T. P., Challinor, V. L., Lehmann, R. P., De Voss, J. J., & Blanchfield, J. T. (2019). Caco-2 monolayer permeability and stability of Chamaelirium luteum (false unicorn) open-chain steroidal saponins. ACS Omega. 4(4), 7658–7666. https://doi.org/10.1021/acsomega. 9b00496 Tulloch, L. B., Martini, V. P., Iulek, J., Huggan, J. K., Lee, J. H., Gibson, C. L., Smith, T. K., Suckling, C. J., & Hunter, W. N. (2010). Structurebased design of pteridine reductase inhibitors targeting African sleeping sickness and the Leishmaniases. Journal of Medicinal Chemistry, 53(1), 221–229. Wang, J., Wang, W., Kollman, P. A., & Case, D. A. (2001). Antechamber: An accessory software package for molecular mechanical calculations. Journal of the American Chemical Society, 222, U403. Wang, J., Wolf, R. M., Caldwell, J. W., Kollman, P. A., & Case, D. A. (2004). Development and testing of a general amber force field. Journal of Computational Chemistry, 25(9), 1157–1174. https://doi.org/10.1002/ jcc.20035 Zikri, A. T., Pranowo, H. D., & Haryadi, W. (2020). Stability, hydrogen bond occupancy analysis and binding free energy calculation from flavonol docked in DAPK1 active site using molecular dynamic simulation approaches. Indonesian Journal of Chemistry, 21(2), 383–390. https://doi.org/10.22146/ijc.56087