Structure-based de novo design, molecular docking and molecular dynamics of primaquine analogues acting as quinone reductase II inhibitors

Accepted Manuscript Title: Structure-Based de Novo Design, Molecular Docking and Molecular Dynamics of Primaquine Analogues Acting as Quinone Reductase II Inhibitors Author: Erika Murce Teobaldo Ricardo Cuya Guizado Helmut Isaac Padilla-Chavarria Tanos Celmar Costa Franc¸a Andre Silva Pimentel PII: DOI: Reference:

S1093-3263(15)30061-9 http://dx.doi.org/doi:10.1016/j.jmgm.2015.10.001 JMG 6610

To appear in:

Journal of Molecular Graphics and Modelling

Received date: Revised date: Accepted date:

8-4-2015 29-9-2015 1-10-2015

Please cite this article as: Erika Murce, Teobaldo Ricardo Cuya Guizado, Helmut Isaac Padilla-Chavarria, Tanos Celmar Costa Franc¸a, Andre Silva Pimentel, StructureBased de Novo Design, Molecular Docking and Molecular Dynamics of Primaquine Analogues Acting as Quinone Reductase II Inhibitors, Journal of Molecular Graphics and Modelling http://dx.doi.org/10.1016/j.jmgm.2015.10.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Structure-Based de Novo Design, Molecular Docking and Molecular Dynamics of Primaquine Analogues Acting as Quinone Reductase II Inhibitors Erika Murcea, Teobaldo Ricardo Cuya Guizadob, Helmut Isaac Padilla-Chavarriaa, Tanos Celmar Costa Françab,c, Andre Silva Pimentela* a

Departamento de Química, Pontifícia Universidade Católica do Rio de Janeiro, RJ 22453-900 Brazil.

b

Laboratory of Molecular Modeling Applied to the Chemical and Biological Defense (LMCBD), Military Institute of Engineering, 22290-270,

Rio de Janeiro, RJ, Brazil. c

Center for Basic and Applied Research, Faculty of Informatics and Management, University of Hradec Králové, Hradec Králové, Czech

Republic. *Corresponding author.

Graphical Abstract

Highlights Five C-6 substituted primaquine new analogues were investigated suggesting that all analogues bind better to the quinone reductase II than primaquine and may become better antimalarials with a reduced toxicity related to the C-5 analogues.

Abstract Primaquine is a traditional antimalarial drug with low parasitic resistance and generally good acceptance at higher doses, which has been used for over 60 years in malaria treatment. However, several limitations related to its hematotoxicity have been reported. It is believed that this toxicity comes from the hydroxylation of the C-5 and C-6 positions of its 8-aminoquinoline ring before binding to the molecular target: the quinone reductase II (NQO2) human protein. In this study we propose primaquine derivatives, with substitution at position C-6 of the 8aminoquinoline ring, planned to have better binding to NQO2, compared to primaquine, but with a reduced toxicity related to the C-5 position being possible to be oxidized. On this sense the proposed analogues were suggested in order to reduce or inhibit hydroxylation and further oxidation to hemotoxic metabolites. Five C-6 substituted primaquine analogues were selected by de novo design and further submitted to docking and molecular dynamics simulations. Our results suggest that all analogues bind better to NQO2 than primaquine and may become better antimalarials. However, the analogues 3 and 4 are predicted to have a better activity/toxicity balance.

Keywords: Primaquine; malaria; quinone reductase II; toxicity; de novo design

Introduction Malaria is an infectious disease caused by the Plasmodium family of protozoa and vectored by mosquitoes. Mainly affecting tropical regions, it is considered an endemic disease in many low-income countries, infecting millions of people and negatively impacting their economic and social development. Despite the successful reduction in the number of cases worldwide, especially due to investments on improving the life quality of African populations and expanding the types of therapies, the disease is still responsible for almost half a million child deaths annually [1]. The development of resistance by the parasites towards the existing treatments urges the discovery of alternative antimalarials [2-4]. Primaquine, (RS)-N-(6-methoxyquinolin-8-yl)pentane-1,4-diamine, is a member of the 8-aminoquinoline group of compounds, commonly used in the treatment and prevention of malaria [5]. First synthesized in 1946, not only it showed less toxicity then the previously used anti-malarial, named as pamaquine, but also acted as a gametocytocide against Plasmodium vivax and Plasmodium ovale, as well as a tissue-schizontocide. It shows efficient prophylactic action, both primarily and terminally, as an anti-relapse treatment. Co-administration with other antimalarials, such as chloroquine [6], maximizes the efficiency of the treatment, especially on the schizontocidal front. This combination is strongly used in endemic areas and towards radical cures. However, parasitic resistance towards chloroquine, a long-standing issue, limits the treatment [7,8]. Low parasitic resistance and generally good acceptance at higher doses justifies why primaquine has been used for over 60 years in malaria treatment [5]. However, several limitations related to its hematotoxicity are present, such as the possibility of inducing hemolytic anaemia, and the increase of methemoglobines in blood flow, which at high doses leads to cyanosis and further death [9]. Additionally, some specific conditions aggravate even more its toxicity, such as the genetic deficiency of glucose-6-phosphate dehydrogenase (G6PD) [9, 10]. Besides, its low plasma half-life urges frequent administration, which amplifies adverse effects. It is known that primaquine is metabolised by cytochrome P450 [11-13] yielding several compounds. Its main metabolite is carboxyprimaquine, a less active and less toxic form of the parent molecule. However, metabolic hydroxylation is responsible for generating hydroxylated species associated to its hematotoxicity due to their redox-cycling which forms hydrogen peroxide and yields free-radicals [14, 15]. The most representative members of this group are the 5-hydroxyprimaquine [5,16-18], 5,6-dihydroxyprimaquine, 6-methoxy-5-hydroxy-8-

aminoquinoline [5], 5,6-dihydroxy-8-aminoquinoline [5], and the 6-methoxy-8-aminoquinoline (MAQ) [5,19-22], which account for most of the drug toxicity. This suggests that the C-5 and C-6 positions of the 8-aminoquinoline ring are the most susceptible ones to suffering hydroxylation. Figure 1 shows the primaquine structure with the ring numbering indicated. The evidence reported in literature [5,14-22] is that the actual active species of primaquine may be the oxidised analogues. These oxidised analogues are probably the same that cause most toxicity. Thus, considering plausible the hypothesis that the oxidised analogues are responsible only of the toxicity and not of the activity of primaquine analogues, it is possible that designed structures with substitutions in only one position (C-5 or C-6) may be oxidised, leaving the rest to the aromatic ring untouched and being active. Anyway, the most important point here is to develop a new prospective primaquine analogue with a better activity/toxicity balance aimed in 6 decades of development in the literature.5 The quinone reductase II (NQO2) human protein is a recently discovered target of primaquine [23]. NQO2 is a homodimer flavindependant reductase, containing a flavin adenine dinucleotide (FAD) tightly bound as a cofactor in each subunit, as well as metal binding sites with zinc atoms. It was shown that primaquine competitively inhibits this protein with respect to the electron donor dihydronicotinamide riboside (NRH), responsible for the reduction of FAD to FADH2. Primaquine is only able to bind to the oxidized and most stable form of the protein, interacting by - stacking with the neutral isoalloxazine ring of FAD. This interaction changes the NQO2 conformation, suggesting a mechanism of flavin redox switch, where the FAD oxidation state is responsible for regulating the protein interactions [24]. The structure of the oxidized NQO2 complexed with primaquine (NQO2ox-PQ) (PDB code: 4FGJ) was obtained by crystallographic data [25]. This study aims to propose molecules with better binding compared to the primaquine binding with NQO2, following the four ADME criteria that influence the compound levels and kinetics of compound exposure to the tissues, and the toxic potential (endocrine and metabolic disruption, some aspects of carcinogenicity and cardiotoxicity) of compounds. As changes in ring positions C-5 and C-6 of the 8-aminoquinoline ring may reduce the primaquine toxicity, [26] the proposed analogues were suggested in order to reduce or inhibit hydroxylation and further oxidation of primaquine to its hydroxylated, hemotoxic metabolites. The NQO2ox-PQ complex was referred as a model in order to study the

interactions between the proposed molecules and the target protein. Thus, we propose new primaquine analogues to be tested as NQO2 inhibitors. Methodology De novo design Substitutions on the 8-aminoquinoline ring [27] were used in this study as a strategy to reduce toxicity. For this approach, de novo drug design was applied using the LigBuilder 2.0 package [28]. This program uses the three-dimensional structure of a target protein to build ligands iteratively from a library of organic fragments, suggesting hundreds of possible leads. A “seed‖ structure of the aminoquinoline ring was designed as a starting point for the program, using the primaquine structure from PDB (PDB code: 1PQ). The cavity module of LigBuilder 2.0 was employed to put the protein coordinates into lattice space, detect the cavity location and boundary and analyse the physicochemical properties of the binding site. Following this procedure, the build algorithm was applied to generate new molecules according to the shape and properties of the NQO2 binding site [29]. In this case, fragments from the database were added to the pre-selected positions, C-5 and C-6, of primaquine, and tested in the corresponding pocket of NQO2 as showed in Figure 2. In the case of the methoxy group at C-6 position, the fragment was exchanged with the methoxy group at this postion. The pre-defined parameters of the program were modified according to Lipinski rule of 5, in order to optimize the drug likeness of the compounds [30]. A total of 100 runs in the build module were performed in order to sample the solution space. For each procedure, the population was set to 10,000 and 10 generations were carried out. The obtained results were organized in clusters by an empirical score function, which considers binding affinity, biological availability, shape complementarity, and synthesizability in order to rank the results [31]. These results were further processed, filtered, recommended and re-clustered by the program in order to select optimal drug candidates, with higher success rates. The yielded compounds were then submitted to docking studies.

Docking studies

The docking energies were obtained using the software Molegro Virtual Docker (MVD) using the MolDock Score algorithm, an adaptation of the Differential Evolution (DE) algorithm [32]. The MolDock Score energy, Escore, is defined by Equation (1), where Einter is the ligand-protein interaction energy and Eintra is the internal energy of the ligand. Einter is calculated according to Equation (2). Escore =Eint er +Eint ra Eint er =

∑

∑

(1)

[ E PLP (r ij )+332. 0

qi q j ] 4rij 2

i=ligant j=protein

(2)

The EPLP term is a ―piecewise linear potential‖ which uses two different parameters, one for the approximation of the steric term (van der Waals) between atoms, and another, for the potential of hydrogen bonds. The second term describes the electrostatic interactions between charged atoms based on a distance 4rij, where rij represents the distances between the considered atoms [32]. Eintra is calculated according to Equation (3). Eint ra =

∑

∑

i=ligant j=ligand

[ E PLP (r ij )]+

∑

flexiblebonds

A [1− cos(mθ− θ o )]+E clash

(3)

The first term in Equation 3 calculates all the energies involving pairs of atoms of the ligand, except those connected by two bonds. The second term refers to the torsional energy, where θ is the torsional angle of the bond and m and A are parameterized according to the hybridization of the bonded atoms. The average of the torsional energy bond contributions is used when several torsions can be determined. The last term, Eclash, assigns a penalty of 1,000 kcal mol-1 if the distance between two heavy atoms (more than two bonds apart) is smaller than 2.0 Å, ignoring infeasible ligand conformations [32]. The energies of Equations (1)–(3) are given in kcal mol-1. The binding sites of NQO2 are predicted with the grid-based cavity program, which employs a van der Waals molecular surface. The docking is constrained by a sphere with radius of 12 Å, centered in the binding site of primaquine, with a resolution of 0.3 Å. Residues within range of the steric contact for the primaquine molecule are flexible, with a tolerance of 0.9 Å, which refers to the optimal energy of interaction of residues between 3.6 Å and 4.5 Å of the ligand. Due to the stochastic nature of the docking algorithm, the protocol is performed with about 30 runs with an initial population of 50 and maximum of 1500 iterations. For analysis of each analogue, 30 poses are returned for further investigation. As ranked by the scoring function (MolDock score), the conformation of each ligand with the best binding energy is selected for the MD. The docking accuracy of the MolDock score has been

evaluated extensively for many protein targets, and this score is able to identify the correct binding mode of most complexes [32]. Besides, the docking protocol used is validated by re-docking studies.

Molecular dynamics The docking results of the lowest energetic conformations of primaquine and its analogues obtained were submitted to steps of molecular dynamics (MD) simulations in order to study the dynamic behavior of the systems. The MD simulations were carried out using the GROMACS 4.5.5 package [33] in boxes of 5.9x5.4x7.1 nm containing around 16,000 water molecules. The system also contained two zinc ions, glycerol molecules and FAD, in order to accurately simulate primaquine and its analogues interactions with the target. Glycerol and zinc ions are far away from the binding site and do not interact with the ligands, only FAD. The parameterizations of ligands, cofactor, FAD, and glycerol were carried out in the Automated Topology Builder (ATB) and Repository version 2.1 [34] using the GROMOS 53a6 force field [35]. The system was neutralized with 8 Na+ ions. Each system was put into a space-filling box, which undergoes replication periodically in the x, y and z directions, in order to be simulated. These systems were studied using the single-point charge water model [36]. The energy minimization algorithms used were steepest descent with position restraint (PR) of the ligands, followed by steepest descent without PR, conjugate gradients, and finally, a quasi-Newton–Raphson method, using the l-BFGS algorithm [37]. The energy minimized complexes were then subjected to MD simulations in two steps using a NpT ensemble with the velocity rescaling scheme to simulate the human body system at 310 K (τT = 0.1 ps), and the Berendsen pressure coupling scheme to simulate a constant and isotropic pressure at 1 atm (τP = 1 ps) [38]. Initially, we carried out 500 ps of MD at 298K with PR for the protein, in order to ensure a balance of the solvent molecules around the residues of the protein. Subsequently, an MD simulation of 20 ns at 310K without any constraint was carried out using 2 fs of integration time and a cut-off of 10 Å for long-distance interactions which were taken into account by means of the Particle Mesh Ewald (PME) technique [39]. A total of 20,000 conformations were obtained during each simulation. Plots of variation of total energy, distance, variation of root-mean-square deviation (RMSD) and hydrogen

bonds formed along the MD simulations were generated. The frames of MD simulations were generated in the VMD program [40]. All analyses were carried out using coordinates recorded at intervals of 2 ps. The binding energies of the compounds docked into NQO2 were estimated by using the molecular mechanics Poisson−Boltzmann surface area method (MM-PBSA) [41]. The auto correlation time of the protein-ligand system was about 450 ps. Thus, the trajectory was sampled at each 500 ps, in order to avoiding the correlation of each frame. MM-PBSA, as implemented by Kumari et al [41], is designed to calculate the binding energies of small and macro-molecules such as proteins, nucleic acids, and other complex systems. This method combines three contributions to the free energy: the first one corresponds to the potential energy in vacuum, and includes bonded terms (such as bond, angle, and torsion energies) as well as non-bonded terms (such as van der Waals and electrostatic interactions). The solvation effects are considered in the second term, which includes the sum of two energy terms, i.e., the polar and non-polar solvation energies, using an implicit solvation model. The last term includes the entropic effect with the complex in gas phase [41]. The protein-ligand interactions involve several amino acids. In order to quantify and determine the relevance of each amino acid in relation to the ligand, interactions like hydrogen bonding and salt bridges were evaluated. Intermolecular surface contact (ISC) was employed to evaluate the relevance of the residues in protein-ligand interaction [42]. The ISC provides information on the interaction between the ligands and given NQO2 residues. The ISC was obtained by means of the SURF software sample over 1,000 MD frames. This software is based on the Connolly algorithm [43], using a 1.4 Å test radius and 3 point/Å. The ISC is determined by the intersection of the surfaces accessible to the solvent of both protein and ligand.

Toxic potential The five analogues and primaquine as reference were submitted to an in silico tool (VirtualToxLab) [44] for predicting their toxic potential (endocrine and metabolic disruption, some aspects of carcinogenicity and cardiotoxicity). The toxic potential (ToxPot) is based solely on thermodynamic considerations and derived from the binding affinities (K) computed by means of 4D Boltzmann scoring of all binding modes

identified through automated, flexible docking to the 3D structure of 16 target protein (nuclear receptors I (AR, ERα, ERβ, GR, MR, PR), nuclear receptors II (LXR, PPARγ, TRα, TRβ), cytochrome P450 enzymes (1A2, 2C9, 2D6, 3A4) and two singular families (AhR, hERG)) known to trigger adverse effects. The ToxPot is calculated through a non-linear function, ranges from none toxicity (0.0) to extreme toxic (1.0), using a normalization process presented as follows: K > 10-2 M

Knorm = 0.0

10-2 M > K > 10-10 M

Knorm=[log(10−2)−log(K)]/[log(10−10)−log(10−2)]

Affinity < 10-10 M

Knorm = 1.0 (4)

Then, the individual toxic potential (TPindividual) is calculated for each individual target protein in Eq. (5): TPindividual = Knorm x weightsd where weightsd is 1.0 – 0.125x(standard deviation/K). Finally, the TPindividual is ranked by their value and the contribution to the overall toxic potential (TPoverall) is summed by using Eq. (6): TPoverall   1  TPoverall, current  TPindividual, n  Wsuperfamily 16

n 1

(6) where Wsuperfamily=1/n (n is nth member of a super family). When the ToxPot is less than 0.1, this indicates that the compound is unlikely to bind to any of the 16 target proteins, however compounds associated with a ToxPot higher than 0.6 may trigger severe adverse effects or even toxic responses.

Results and Discussion

De novo design

Around 1,000 compounds were obtained and grouped together in two clusters ordered by the score function, with average similarity of 96% and 95%, respectively. These results were submitted to processing, filtering, score calculation, recommendation and re-clustering by the build algorithm. After the processes, two compounds were selected from the first cluster (analogues 1 and 2), and three from the second cluster (analogues 3, 4 and 5) with the following properties of the Lipinski rule of 5, as presented in Table 1. Their chemical structures are presented in Figure 3. All five compounds featured fragments added to the position C-6 of the 8-aminoquinoline ring connected to the main ring by amide and carbamide fragments. The selected compounds featured amine (NH2), hydroxyl (OH), and hydroxymethyl (CH2OH) groups in orto, meta and para positions. Despite the fact that primaquine analogues with fragments substituted in the C-5 position showed lower toxicity than primaquine in the literature [45], our study did not find any compound with fragments in this position.

Docking

Figure 4 shows the re-docking of primaquine on its experimental structure (PDB: 4FGJ). The RMSD generated by the docking protocol used in this work was 1.27 Å. Keeping in mind that a RMSD value under 2.000 Å is considered acceptable [46-48], these results validate the docking protocol used. After that the five analogues and primaquine were, then, submitted to flexible molecular docking inside NQO2 in order to evaluate the differences between the analogues and the parent drug in terms of binding energy to the target. It was observed that all five analogue-target complexes presented higher stability than the primaquine-macromolecule system. Results in Table 2 show that analogue 5 presents both the best binding and H-bond energies contributions. The ligands conformations with the lowest energy values were selected for the MD simulations.

Molecular dynamics

The dynamic behaviours of the systems were studied considering the RMSD of each analogue towards the protein throughout the 20 ns of MD simulation. The RMSD calculation evaluates the stability of the compounds inside the target cavity over time in dynamic conditions. As can be seen in Figure 5, all analogue-protein complexes showed lower fluctuation in RMSD than the complex primaquine-protein. This result suggests stronger interactions between the analogues and the target when compared to primaquine.

Table 3 shows the prevalence of H-bonds obtained from the MD data. It can be seen that analogues 2, 4 and 5 show the highest values of prevalence, corroborating the H-bond energies observed in the docking studies (see table 2). The analogue-protein complexes were stabilised by several hydrogen bond interactions along the MD simulations. Approximately two times more interactions were observed for the analoguesprotein complexes than for the primaquine-protein complex. This behaviour was already expected due to the presence of CH2OH, NH2 and OH groups in the phenyl fragments added by the de novo design. Our results also show that Asn161 act as donor and acceptor for analogue 2 while Trp105 is a donor and Asp117 is an acceptor for analogue 4. For analogue 5, Gln122 and Trp132 are donors and Asp117 is an acceptor. These results suggest that analogue 5 has a more extended and strong H-bond network compared with the other analogues. It was also observed that all analogues make H-bonds with the cofactor FAD, during the MD simulations but with less than 1.0% of prevalence. However, - stacking interactions were continuously observed between the analogues and cofactor, but analogues do not interact with the zinc atom.

Figure 6 show the plots of ISC for the analogues inside NQO2. Here only amino acids with ISC > 10 were considered relevant for the study [42]. It can be observed that Phe126, Glu193 and Gln122 are important residues in the ligands-target interactions. However, some residues such as Phe, Met, Trp, Gln and Tyr around the active site seem to be interacting with all analogues but not with primaquine, supporting our previous results suggesting that they bind more tightly inside the pocket than primaquine. Thus, these interactions show an important role in stabilizing the analogues-protein complexes.

The MM-PBSA results, reported in Figure 7, corroborate with the analysis performed above and also point to the five analogues as stronger binders inside NQO2 than primaquine. The binding energy of analogue 5 with the NQO2 binding site is approximately three times larger than that found for primaquine, making the analogue 5-NQO2 complex more stable. Considering that MM-PBSA is a relative binding energy calculation method, we compared the results to the binding energies of the ligands obtained using the MVD [32]. Results, shown in Figure 8, suggest a good correlation between both energy methods with a R2 = 0.9.

Toxic potential Table 4 presents the logarithm of binding affinities of the five analogues and primaquine as reference to sixteen target protein. Also, it shows the toxic potential of each analogue and primaquine. The analogue 1 belongs to the toxic potential class III, which strongly bind to two target classes. The analogues 2 and 5 moderately bind to one target class, belonging to the toxic class I. Finally, the compounds 3 and 4 belong to the same toxic potential class of primaquine, class 0, whose compounds modestly bind to a single target classor weakly to several target classes. These latter analogues have lower binding to the cytochrome P450 enzymes (1A2, 2C9, 2D6, 3A4) compared to primaquine, therefore may have a better activity/toxicity balance aimed in the literature.5

Conclusion De novo drug design was successfully applied to propose new primaquine analogues. Then, these analogues were submitted to molecular docking and dynamics simulations to understand the interactions of analogues with the NQO2 human protein. Our results suggest that the five primaquine analogues with substitutions on position C-6 of the 8-aminoquinoline ring, selected by de novo design, interact better with NQO2 than primaquine. Considering that hydroxylation on this position is directly related to toxicity of primaquine, it is reasonable that these analogues may be less toxic. Our docking and MD achievements are in agreement and the MM-PBSA results are highly correlated with the molecular docking energies. It was found that Phe126, Glu193 and Gln122 strongly interact with the analogues, however amino acids Phe, Met, Trp, Gln and Tyr around the active site also contribute significantly to the analogue-target binding. The analysis performed pointed to analogue 5 as the

one with the most stable interactions inside NQO2 binding site, among the molecules studied. The analogues 3 and 4 have promising predictions on their toxicity compared with primaquine, intermediate binding interaction, but better than primaquine. They seem to have a weaker interaction with the cytochrome P450 enzymes, which may produce a better activity/toxicity balance compared with primaquine. Thus, this molecule seems to be a promising candidate to be tested as a potential new antimalarial drug.

Acknowledgments

The authors thank the Brazilian funding agencies CNPq, CAPES (Grant No. 02559/09-9), and FAPERJ (Grant E-26/102.993/2012) for financial support. T.R.C.G. also thanks CAPES for a post-doc PNPD fellowship. H. I. P. C., E.M.S. and T.C.C.F. acknowledge CNPq for research fellowships (Grants 304557/2012-9 and 400856/2014-0). This work was also supported by Excellence project FIM.

References [1] World Health Organization (2014). World malaria report 2014. [2] H. Ocholla, M. D. Preston, M. Mipando, A. T. Jensen, S. Campino, B. MacInnis, D. Alcock, A. Terlouw, I. Zongo, J. B. Oudraogo, A. A. Djimde, S. Assefa, O. K. Doumbo, S. Borrmann, A. Nzila, K. Marsh, R. M. Fairhurst, F. Nosten, T. J. Anderson, D. P. Kwiatkowski, A. Craig, T. G. Clark, Montgomery J. Whole-genome scans provide evidence of adaptive evolution in Malawian Plasmodium falciparum isolates. J Infect Dis, 210(12) (2014), pp. 1991-2000. doi: 10.1093/infdis/jiu349. [3] R. N. Price, L. von Seidlein, N. Valecha, F. Nosten, J. K. Baird, N. J. White. Global extent of chloroquine-resistant Plasmodium vivax: a systematic review and meta-analysis. Lancet Infect Dis, 14(10) (2014), pp. 982-991. doi: 10.1016/S1473-3099(14)70855-2. [4] S. Sinha, B. Medhi, R. Sehgal. Challenges of drug-resistant malaria. Parasite, 21 (2014), pp. 61. doi: 10.1051/parasite/2014059. [5] N. Vale, R. Moreira, P. Gomes. Primaquine revisited six decades after its discovery. Eur J Med Chem, 44(3) (2009), pp.937-953. doi: 10.1016/j.ejmech. [6] J. Ducharme, R. Farinotti. Clinical pharmacokinetics and metabolism of chloroquine. Focus on Recent Advancements. Clin. Pharmacokinet, 31 (1996), pp. 257–274. doi: 10.2165/00003088-199631040-00003. [7] D. C. Warhurst, R. Killick. Spontaneous resistance to chloroquine in a strain of rodent malaria (plasmodium berghei yoelii). Nature, 213(5080) (1967), pp 1048-&. doi:10.1038/2131048a0. [8] C. D. Fitch. Chloroquine resistance in malaria: a deficiency of chloroquine binding. Proc Natl Acad Sci USA, 64(4) (1969), pp. 1181-7. [9] T. S. Nascimento, R. Pereira, H. Mello, J. Costa. Metemoglobinemia: do Diagnóstico ao Tratamento. Rev Bras Anestesiol, 58(6) (2008), pp. 651-664. doi: 10.1590/S0034-70942008000600011. [10] P. Kedar, P. Warang, S. Sanyal, R. Devendra, K. Ghosh, R. Colah. Primaquine-induced severe methemoglobinemia developed during treatment of Plasmodium vivax malarial infection in an Indian family associated with a novel mutation (p.Agr57Trp) in the CYB5R3 gene. Clin Chim Acta, 437 (2014), pp. 103-105. doi: 10.1016/j.cca.2014.07.015.

[11] B. S. Pybus, J. C. Sousa, X. Jin, J. A. Ferguson, R. E. Christian, R. Barnhart, C. Vuong, R. J. Sciotti, G. A. Reichard, M. P. Kozar, L. A. Walker, C. Ohrt, V. Melendez. CYP450 phenotyping and accurate mass identification of metabolites of the 8-aminoquinoline, anti-malarial drug primaquine. Malaria J, 11 (2012), 259. doi: 10.1186/1475-2875-11-259. [12] B. S. Pybus, S. R. Marcsisin, X. Jin, G. Deye, J. C. Sousa, Q. Li, D. Caridha, Q. Zeng, G. A. Reichard, C. Ockenhouse, J. Bennett, L. A. Walker, C. Ohrt, V. Melendez. The metabolism of primaquine to its active metabolite is dependent on CYP 2D6. Malaria J, 12(212) (2013). doi: 10.1186/1475-2875-12-212. [13] S. Ganesan, B. L. Tekwani, R. Sahu, L. M. Tripathi, L. A. Walker. Cytochrome P-450-dependent toxic effects of primaquine on human erythrocytes. Toxicol Appl Pharm, 241 (2009), pp. 14-22. doi: 10.1016/j.taap.2009.07.012. [14] J. Vásquez-Vivar, O. Augusto. Hydroxylated metabolites of the antimalarial drug primaquine - Oxidation and redox cycling. J Biol Chem, 267(10) (1992), pp. 6848-6854. [15] J. Vásquez-Vivar, O. Augusto. Oxidative activity of primaquine metabolites on rat erythrocytes IN vitro and in vivo. Biochem pharmacol, 47(2) (1994), pp 309-316. [16] Z. S. Bowman, D. J. Jollow, D. C. McMillan. Primaquine-induced hemolytic anemia: Role of splenic macrophages in the fate of 5hydroxyprimaquine-treated rat erythrocytes. J Pharmacol Exper Ther, 315(3) (2005), pp. 980-986. doi: 10.1124/jpet.105.090407. [17] Z. S. Bowman, J. D. Morrow, D. J. Jollow, D. C. McMillan. Primaquine-induced hemolytic anemia: Role of membrane lipid peroxidation and cytoskeletal protein alterations in the hemotoxicity of 5-hydroxyprimaquine. J Pharmacol Exper Ther, 314(2) (2005), pp. 838-845. doi: 10.1124/jpet.105.086488. [18] Z. S. Bowman, J. E. Oatis Jr, J. L. Whelan, D. J. Jollow, D. C. McMillan. Primaquine-induced hemolytic anemia: Susceptibility of normal versus glutathione-depleted rat erythrocytes to 5-hydroxyprimaquine. J Pharmacol Exper Ther, 309(1) (2004), pp.79-85. doi: 10.1124/jpet.103.062984. [19] J. D. Baty, D. A. Evans, P. A. Robinson. Identification of 6-methoxy 8-aminoquinoline as a metabolite of primaquine in man. Biom Mass Spectrom, 2(6) (1975), pp. 304-306. doi: 10.1002/bms.1200020605.

[20] L. J. Bolchoz, A. K. Gelasco, D. J. Jollow, D. C. McMillan. Primaquine-induced hemolytic anemia: Formation of free radicals in rat erythrocytes exposed to 6-methoxy-8-hydroxylaminoquinoline. J Pharmacol Exper Ther, 303(3) (2002), pp. 1121-1129. doi: 10.1124/jpet.102.041459. [21] L. J. Bolchoz, J. D. Morrow, D. J. Jollow, D. C. McMillan. Primaquine-induced hemolytic anemia: Effect of 6-methoxy-8hydroxylaminoquinoline on rat erythrocyte sulfhydryl status, membrane lipids, cytoskeletal proteins, and morphology. J Pharmacol Exper Ther, 303(1) (2002), pp. 141-148. doi:10.1124/jpet.102.036921. [22] L. J. Bolchoz, R. A. Budinsky, D. C. McMillan, D. J. Jollow. Primaquine-induced hemolytic anemia: formation and hemotoxicity of the arylhydroxylamine metabolite 6-methoxy-8-hydroxylaminoquinoline. J Pharmacol Exper Ther, 297(2) (2001), pp. 509-515. [23] P. R. Graves, J. J. Kwiek, P. Fadden, R. Ray, K. Hardeman, A. M. Coley, M. Foley, T. A. Haystead. Discovery of novel targets of quinoline drugs in the human purine binding proteome. Mol Pharmacol, 62(6) (2002), pp. 1364-1372. doi: 10.1124/mol.62.6.1364. [24] J. J. Kwiek, T. A. Haystead, J. Rudolph. Kinetic mechanism of quinone oxidoreductase 2 and its inhibition by the anitmalarial quinolines. Biochemistry-US, 43(15) (2004), pp. 4538-4547. doi: 10.1021/bi035923w. [25] K. K. Leung, B. H. Shilton. Chloroquine Binding Reveals Flavin Redox Switch Functin of Quinone Reductase 2. J Biol Chem, 288 (2013), pp. 11242-11251. doi: 10.1074/jbc.M113.457002. [26] C. Teixeira, N. Vale, B. Pérez, A. Gomes, J. Gomes, P. Gomes. ―Recycling‖ classical drugs for malaria. Chem Rev, 114 (2014), pp. 1116411220. doi 10.1021/cr500123g. [27] R. G. Ridley, A. T. Hudson. Quinoline Antimalarials. Expert Opin Ther Pat, 8(2) (1998), pp. 121-136. doi: 10.1517/13543776.11.2.185 [28] Y. Yuan, J. Pei, L. Lai. LigBuilder 2: A Practical de Novo Drug Design Approach. J Chem Inf and Model, 51 (2011), pp. 1083–1091. doi: 10.1021/ci100350u. [29] S. Ni, Y. Yuan, J. Huang, X. Mao, M. Lv, J. Zhu, X. Shen, J. Pei, L. Lai, H. Jiang, J. Li. Discovering Potent Small Molecule Inhibitors of Cyclophilin A Using de Novo Drug Design Approach. J Med Chem, 52 (2009), pp. 5295–5298. doi: 10.1021/jm9008295.

[30] C. A. Lipinski, F. Lombardo, B. W. Dominy, P. J. Feeney. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv Drug Deliver Rev, 23 (1997), pp. 3-25. doi:10.1016/S0169-409X(96)00423-1. [31] R. X. Wang, L. Liu, L. Lai, Y. Tang. SCORE: A new empirical method for estimating the binding affinity of a protein-ligand complex. J Mol Model, 4(12) (1998), pp. 379-394. doi: 10.1007/s008940050096. [32] R. Thomsen, M. H. Christensen. MolDock: a new technique forhigh-accuracy molecular docking, J Med Chem 49 (2006), pp. 3315–3321. doi:10.1021/jm051197e. [33] B. Hess, C. Kutzner, D. van der Spoel, E. Lindhal. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J Chem Theory Comp 4(3) (2008), pp. 435–447. doi: 10.1021/ct700301q. [34] K. B. Koziara, M. Stroet, A. K. Malde, A. E. Mark. Testing and validation of the Automated Topology Builder (ATB) version 2.0: Prediction of hydration free enthalpies. J Comput Aided Mol Des, 28(3) (2013), pp. 221-33. doi:10.1007/s10822-014-9713-7. [35] C. Oostenbrink, A. Villa, A. E. Mark, W. F. van Gunsteren. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J Comput Chem, 25(13) (2004), pp. 1656-1676. doi: 10.1002/jcc.20090. [36] H. Weinstein, S. Topiol, R. Osman. On the relation between charge redistribution and intermolecular forces in models for molecular interactions in biology. Struct Bond, (1981), pp. 383–396. [37] C. Zhu, R. H. Byrd, J. Nocedal. Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization, ACM. Trans Math Softw 23 (1997), pp. 550 –560. doi:10.1145/279232.279236. [38] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. Dinola, J. R. Haak. Molecular-dynamics with coupling to an external bath. J Chem Phys, 81 (1984), pp. 3684–3690. doi: 10.1063/1.448118. [39] T. Darden, L. Perera, L. Li, L. Pedersen. New tricks for modelers from the crystallography toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations. Structure 7 (1999), pp. R55–R60. [40] W. Humphrey, A. Dalke, K. Schulten. VMD: Visual molecular dynamics, J Mol Graph 14 (1996), pp. 33 –38. doi:10.1016/02637855(96)00018-5.

[41] R. Kumari, R. Kumar, A. Lynn. g_mmpbsa - A GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inf Model 54 (2014), pp. 1951-1962. doi: 10.1021/ci500020m. [42] T. R. C. Guizado, S. R. Louro, C. Anteneodo. Dynamics of heme complexed with human serum albumin: a theoretical approach. Eur Biophys J. 41 (2012), pp. 1033–1042. doi: 10.1007/s00249-012-0860-2. [43] M. L. Connolly. Solvent-accessible surfaces of proteins and nucleic acids. Science, 221 (1983), pp.709–713. [44] A. Vedani, M. Dobler, Z. Hu, M. Smieško. OpenVirtualToxLab—A platform for generating and exchanging in silico toxicity data. Toxycol Lett 2015, 232, 519-532. [45] H. Liu, B. L. Tekwani, N. P. Nanayakkara, L. A. Walker, R. J. Doerksen. Methemoglobin Generation by 8-Aminoquinolines: Effect of Substitution at 5-Position of Primaquine. Chem Res Toxicol, 26 (2013), pp. 1801-1809. doi: 10.1021/tx400067a. [46] G. L. Warren, A. C. Webster, A. M. Capelli, B. Clarke, J. Lalonde, M. H. Lambert, M. Lindvall, N. Nevins, S. F. Semus, S. Senger, G. Tedesco, I. D. Wall, J. M. Woolven, C. E. Peishoff, M. S. Head. A Critical Assessment of Docking Programs and Scoring Functions. J Med Chem 49 (2006), pp. 5912-31. doi: 10.1021/jm050362n. [47] A. R. Leach, B. K. Shoichet, C. E. Peishoff Prediction of Protein-Ligand Interactions. Docking and Scoring: Successes and Gaps. J Med Chem 49 (2006), pp. 5851-55. doi: 10.1021/jm060999m. [48] M. Kontoyanni, L. M. Mcclellan, G. S. Sokol. Evaluation of Docking Performance: Comparative Data on Docking Algorithms. J Med Chem 47 (2004), 558-65. doi: 10.1021/jm0302997.

Figure Captions Figure 1. The primaquine structure with its ring numbers indicated. Figure 2. Quinone reductase II with the analogue 1 in the catalytic site. The analogue 1 and the cofactor are drawn in licorice and van der Waals styles respectively. Figure 3. The five selected molecules for further studies, as outputted from LigBuilder. Figure 4. Redocking calculation of the primaquine. In yellow the primaquine from the PDB structure 4FGJ, in red the primaquine conformation obtained from the re-docking with a RMSD = 1.27 A. Figure 5. The root mean square deviation (RMSD) of primaquine and all five analogues during the 20 ns molecular dynamics. Figure 6. Intermolecular surface contact of primaquine and all five analogues. Figure 7. MM-PBSA energies of primaquine and all five analogues calculated during the 20 ns of MD simulation. Figure 8. The correlation between MM-PBSA and docking energies for primaquine and the five analogues.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Tables Analogue 1

Analogue 2

Analogue 3

Analogue 4

Analogue 5

Score

6.73

6.47

6.78

6.71

6.43

LigScore

4.13

3.84

4.05

3.98

3.62

Stability

61%

58%

60%

58%

56%

Strong HB

1.98

1.43

1.00

0.95

1.00

Hydrophobicity

2.95

2.95

2.99

2.99

2.66

XLogP

1.87

1.47

1.36

1.05

1.36

Weight

380

394

395

395

395

Table 1. Results for score and Ligscore from Ligbuilder, and the following Lipinski rule of 5 properties for each analogue.

Analogues

MolDock

H-bonds

H-bond energy

Binding energy

Score (kcal/mol)

interactions

analogue 1

-183.47

Glu-193, Gly-68

-7.499

-191.722

analogue 2

-201.068

FAD, Gly-174, Thr-71

-8.288

-225.149

analogue 3

-200.694

FAD, Asp-117

-7.092

-223.23

analogue 4

-199.881

FAD, Thr-71, Asn-161

-8.166

-228.258

analogue 5

-228.441

FAD, Leu-120, Asp-117

-10.368

-241.816

-----

-157.024

(kcal/mol)

Trp-105 primaquine

-149.379

FAD, Gly-68

Table 2. Score, hydrogen bond interactions, H-bond and binding energies obtained from docking calculations.

Analogue 1 Donor

Aceptor

Analogue 2 Prev. (%)

Donor

Aceptor

Analogue 3 Prev. (%)

Donor

Aceptor

Prev. (%)

ASN161ND2

A1-237O1

34.57

ASN161ND2

A2-237O32

97

GLY150N

A3-237O33

10.61

TYR155OH

A1-237O1

35.8

TYR132OH

A2-238O31

7.99

GLN122N

A3-238O31

12.61

ASN161ND2

A1-238O1

30.86

TYR155OH

A2-238O32

15.36

THR71OG1

A3-238O31

19.85

TYR155OH

A1-238O1

44.44

A2-237N33

ASN161OD1

54.06

VAL69N

A3-237O31

16.73

A3-237O31

GLN122OE1

6.49

A3-237N311

GLN122OE1

7.37

A3-238O31

GLY68O

6.37

A3-238O31

ASP117O

11.49

A3-238N311

GLN122OE1

15.11

Analogue 4 Donor

Aceptor

Analogue 5 Prev. (%)

Donor

Aceptor

Primaquine Prev. (%)

Donor

Aceptor

Prev. (%)

TRP105N

A4-238O33

96

PHE126N

A4-238O33

15.58

ASN161ND2

PQ238O1

13.11

TRP105NE1

A4-237O31

14.29

FAD235O10

PQ237N1

48.56

THR71OG1

A4-238O31

42.86

FAD236O10

PQ238N1

8.49

A4-237O31

ASP117O

58.44

PQ237N1

GLU193OE1

18.85

A4-237N21

GLY149O

20.78

PQ237N1

GLU193OE2

14.98

A4-237N311

THR71OG1

24.68

PQ237N1

FAD235O5

13.98

A4-237N311

ASP117O

12.99

PQ238N1

ASN66OD1

8.86

A4-238N34

LEU120O

42.86

PQ238N1

GLN122OE1

9.49

A4-238N311

ASP117O

31.17

Table 3: H-bonds prevalences obtained from MD data. A1, A2, A3, A4 and A5 are referred to the analogues, FAD is referred to the cofactor and PQ is referred to the primaquine. Because the NQO2 is a dimmer, two molecules of each analogue, primaquine and cofactor are show in the table.

Molecule

ToxPot

AR

AhR

1A2

2C9

2D6

3A4

ERα

ERβ

GR

hERG

LXR

MR

PPARОі

PR

TRα

TRβ

Analogue 1

0.750

-6.37

-6.35

*

-5.86

-7.69

-5.94

-7.61

-5.88

-5.91

-6.06

*

-9.76

-5.16

-6.07

*

*

Analogue 2

0.597

-5.40

-8.05

*

-7.05

-7.90

-6.44

-7.85

-4.01

-6.85

-7.41

-4.38

-6.60

-5.61

-4.66

*

-4.97

Analogue 3

0.441

*

-6.01

-5.33

*

-5.56

*

-6.58

*

*

-5.59

*

-5.85

-5.02

*

*

*

Analogue 4

0.458

-6.31

*

-5.13

*

*

*

*

-6.89

*

-4.82

*

-4.76

-4.02

*

*

*

Analogue 5

0.583

-4.95

*

-5.86

-5.26

-7.59

-4.16

-4.48

-8.08

-5.30

-6.12

-4.57

-6.94

-4.77

*

*

*

Primaquine

0.448

-5.28

-6.08

*

-5.08

-5.85

-4.39

-5.49

-6.53

-5.95

-5.84

-4.07

-5.02

-4.47

-5.39

*

*

Table 4. The logarithm of binding affinities computed by means of 4D Boltzmann scoring of all binding modes identified through automated, flexible docking to the 3D structure of the target protein: nuclear receptor class I (AR, ERα, ERβ, GR, MR, PR), nuclear receptor class II (LXR, PPARγ, TRα, TRβ), cytochromes (1A2, 2C9, 2D6, 3A4), hERG, and AhR. Toxpot is a non-linear function of binding affinity ranging from 0.0 (none) to 1.0 (extreme). * means no binding. Dark colors represent high affinities to the target protein, and clear colors the opposity.