- Email: [email protected]
Understanding the binding of daunorubicin and doxorubicin to NADPH-dependent cytosolic reductases by computational methods
European Journal of Medicinal Chemistry 56 (2012) 145e154
Contents lists available at SciVerse ScienceDirect
European Journal of Medicinal Chemistry journal homepage: http://www.elsevier.com/locate/ejmech
Original article
Understanding the binding of daunorubicin and doxorubicin to NADPH-dependent cytosolic reductases by computational methods Davide Pirolli a, Bruno Giardina a, Alvaro Mordente a, Silvana Ficarra b, Maria Cristina De Rosa c, * a
Istituto di Biochimica e Biochimica Clinica, Università Cattolica del Sacro Cuore, Largo F. Vito 1, 00168 Rome, Italy Dipartimento di Chimica Organica e Biologica, Università di Messina, V. le Ferdinando Stagno d’Alcontres 31, 98166 Messina, Italy c CNR-Istituto di Chimica del Riconoscimento Molecolare, c/o Istituto di Biochimica e Biochimica Clinica, Università Cattolica del Sacro Cuore, Largo F. Vito 1, I-00168 Rome, Italy b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 May 2012 Received in revised form 2 August 2012 Accepted 16 August 2012 Available online 24 August 2012
The anthracycline anticancer agents daunorubicin (DAUN) and doxorubicin (DOX) are reduced by different NADPH-dependent cytosolic reductases into their corresponding alcohol metabolites daunorubicinol (DAUNol) and doxorubicinol (DOXol), which have been implicated in the development of chronic cardiomyopathy. To better understand the individual importance of each enzyme in the reduction and to provide deeper insight into the binding at atomic level we performed molecular docking and dynamics simulations of DAUN and DOX into the active sites of human carbonyl reductase 1 (CBR1) and human aldehyde reductase (AKR1A1). Such simulations evidenced a different behavior between the reductases with respect to DAUN and DOX suggesting major contribution of CBR1 in the reduction. The results are in agreement with available experimental data and for each enzyme and anthracycline pair provided the identification of key residues involved in the interactions. The structural models that we have derived could serve as a useful tool for structure-guided drug design studies. Ó 2012 Elsevier Masson SAS. All rights reserved.
Keywords: Anthracycline Cardiotoxicity Molecular docking Molecular dynamics
1. Introduction The anthracycline antibiotics doxorubicin (DOX) and its congener, daunorubicin (DAUN) are among the most effective chemotherapeutic agents in the treatment of numerous malignancies [1]. Since their first clinical trials these drugs have been extensively studied alone and in combination, and they now have a major role in the effective treatment of acute leukemia, nonHodgkin lymphomas, breast cancer, Hodgkin’s disease, and sarcomas [2]. However, two major problems limit their clinical application: drug resistance [3] and cardiotoxicity [4]. This last complication represents a major obstacle to prolonged use of the drug and/or cumulative dose exceeding 550 mg/m2 body surface area [5]. A number of mechanisms have been proposed to account for the anthracycline-induced cardiotoxicity [6]. Among them, a prevailing hypothesis implicates the C-13 hydroxy metabolite, doxorubicinol (DOXol) and daunorubicinol (DAUNol), compounds that lack the antiproliferative activity of DOX and DAUN, as a major component in the development of chronic cardiomyopathy [7]. The
* Corresponding author. E-mail addresses: [email protected] (D. Pirolli), bgiardina@ rm.unicatt.it (B. Giardina), [email protected] (A. Mordente), sficarra@ isengard.unime.it (S. Ficarra), [email protected] (M.C. De Rosa). 0223-5234/$ e see front matter Ó 2012 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ejmech.2012.08.023
accumulation of the metabolite in the heart has been described suggesting that the C-13 hydroxy metabolite might be transported and selectively trapped by cardiac cells [8]. Several cytosolic reductases of the short chain dehydrogenase class (SDR) and of the aldoeketo reductase class (AKR) metabolize DOX to DOXol and DAUN to DAUNol (Fig. 1) implicating these enzymes in the development of cardiotoxicity [9,10]. Enzymes that have actually been described to reduce anthracyclines to C13-hydroxymetabolites are carbonyl reductase 1 (CBR1) [11e15], carbonyl reductase 3 (CBR3) [16], aldo-keto reductase (AKR) 1A1 [12,14,17], AKR1B1 [14,18,19], AKR1B10 [14,18,19], AKR1C3 [14] and AKR1C4 [14]. Divergent results are reported on the ability of AKR1C2 to metabolize doxorubicin and daunorubicin in vitro [12,17]. At present, CBR1 and AKR1A1 are considered the major anthracycline-metabolizing enzymes, at least in human heart [20] and indications about DOX and DAUN reduction obtained from experimental studies are reported in Table 1. Much effort has been devoted toward developing novel anthracyclines showing increased cardiac tolerability with little success. A decrease in the level of alcohol metabolite formation was observed with epirubicin (EPI), the C-40 epimer of DOX, and sabarubicin (MEN 10755), a DOX analog lacking the methoxy group at C-4 and characterized with a disaccharide in the place of daunosamine in C-7 [21]. However these anthracyclines are still characterized by a low but significant cardiotoxicity that makes necessary further investigations toward a less toxic anticancer
146
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 1. Enzymatic reduction of anthracyclines to C-13 secondary alcohol metabolites. Schematic representation of two-electron reduction of the side chain carbonyl group of DOX and DAUN resulting in DOXol and DAUNol, respectively.
anthracycline [22,23]. The evaluation of efficient inhibitors of anthracycline reductases has been also suggested as one potential strategy to optimize cardioprotection [24,25], but structure-based design of specific inhibitors is lacking. On the whole, despite the large amount of available biochemical data, there is no structural information at the atomic level regarding these interactions, how cytosolic reductases bind the anthracyclines within the protein cavity and which residues are involved in the binding process. In this study we have used computational methods, based on molecular docking studies and molecular dynamics (MD) simulation, to explore these features in some detail and to provide useful insights into the future rational design of specific inhibitors of anthracycline reductases activity as well as of novel, less cardiotoxic anthracyclines. Our computational studies correlated with previous knowledge in the field and provide a plausible binding model of DAUN and DOX within the active site of CBR1 and AKR1A1. 2. Results and discussion
2.2. Docking calculations
In carrying out the analysis of DOX and DAUN binding in cytosolic reductases, substrates were first docked into the enzyme active site. In this step the substrate was flexible while the enzyme was fixed in geometry. To account for protein flexibility and to allow the geometry to relax, in a second step the low energy docked conformations were submitted to molecular dynamics simulations. Starting coordinates for CBR1 were obtained from the Protein Data Bank while for AKR1A1 atomic coordinates were generated by homology modeling. 2.1. Construction of AKR1A1 model The structure of human aldehyde reductase at the resolution of 2.48 Å (PDB ID: 2ALR) does not completely resolve the polypeptide chain. In fact, residues 214e225 belonging to the extremely flexible Table 1 Experimental measurements, obtained from the literature, for the reaction of CBR1 and AKR1A1 with the substrates DOX and DAUN. Kinetic parameter
kcat(s1) Km(mM) Vmax/Km (U/mg/M) kcat(min1) Km(mM) kcat(s1) Km(mM) Specific activitya a
CBR1
AKR1A1
Ref
DAUN
DOX
DAUN
DOX
1.27 0.03 108 7 2140 e e 1.94 51 13 e
0.27 0.02 90 16 e 0.62 167 0.21 287 57 e
e e 312 e e e e 120 3
e e e 0.04 247 e e 11 1
nmol/min per mg of protein.
loop B of aldehyde reductase are missing in the 2ALR structure. A 3D structure model of AKR1A was then obtained using homology method. The appropriate template was chosen based on sequence similarity, residue completeness, and crystal resolution. The target and the template, porcine aldehyde reductase (PDB ID: 3H4G), share 94% of their sequences which guarantees the accuracy of the prediction. Among the conformations generated by MODELLER one was selected on the basis of the lowest MODELLER objective function. The Ca root-mean square deviation (RMSD) of the homology model to the template structure was 0.24 Å. The homology model was also compared to the human aldehyde reductase crystal structure (2ALR) comparing only the residues present in the crystal structure. RMSD values were also calculated for heavy atoms within 3.5 Å of the co-crystallized fidarestat (residues 21, 48, 49, 112, 113, 298, 299, 300). The Ca RMSD for the model and the AKR1A crystal structure (2ALR) was 0.9 Å whereas the RMSD of fidarestat binding site residues, within 3.5 Å of the ligand was 1.2 Å. The stereochemical quality of the selected model was examined using the program PROCHECK. The distribution of the F and J angles in the Ramachandran plot is well within the allowed regions. The majority (94%) of the residues are in the most favored regions of the plot; 6% are in the additional allowed regions; no residues lie in disallowed regions. The final structure was further checked by Verify3D, as implemented in the Verify3D Structure Evaluation Server (http://nihserver.mbi.ucla.edu/Verify_ 3D/), a program that analyzes the compatibility of an atomic model (3D) with its own amino acid sequence (1D). As shown in Fig. 2, the VERIFY-3D scores of our model are always positive and similar to those of the template protein and those of human aldehyde reductase. Taken together, these data suggest that built 3D model of AKR1A1 was of reasonable quality compared to 3H4G and 2ALR crystal structures.
[13] [13] [12] [14] [14] [15] [15] [18]
Cytosolic reductases were shown to follow a sequential ordered mechanism in which NADPH binds to the enzyme first and NADPþ is released last after the alcohol product. The reaction mechanism involves a hydride transfer from NADPH to the carbonyl carbon of the substrate and a proton transfer from a Tyr residue at the bottom of the active site cavity, which serves as the general acid catalyst, to the forming alcohol group. Additional active site residues form a proton relay system and stabilize the orientation of the substrate. The conformations from the docking experiments were therefore analyzed for a catalytically productive orientation of the substrate within the binding site. Two criteria were used to identify productive conformations: 1) the distance between the C13 anthracycline carbon atom and the C4 atom of NADPH should be short enough to allow efficient hydride transfer; 2) the oxygen of C13 carbonyl should be located in proximity to the Tyr residue which acts as proton donor during the catalysis. 2.2.1. Test of the method The structures of CBR1 bound to hydroxy-PP and porcine aldehyde reductase bound to fidarestat were chosen to test the docking protocol. The results of docking validation indicate that the predicted binding conformations determined by AutoDock Vina match well with those of the co-crystallized ligand. 2.2.2. Binding modes of DAUN and DOX with CBR1 The overall structure of CBR1 retains the characteristic a/ b folding pattern with a Rossmann-fold of the SDR family and three residues (Ser139, Tyr193 and Lys197) were proposed as the catalytic triad in this short-chain dehydrogenase/reductase enzyme [26,27]. Accordingly, Tyr193 functions as the catalytic acid/base,
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
147
Fig. 2. Average Verify3D score profile. Plots obtained by the Verify-3D program, with a sliding 21-residue window, using the coordinates of AKR1A1 homology model (black), the template protein 3H4G (black dashed) and the crystallographic aldehyde reductase 2ALR (gray). Note that residues 214e224 were unresolved in the 2ALR structure and are therefore missing in the corresponding plot.
Ser139 stabilizes the substrate by forming interactions to the substrate carbonyl, and Lys197 forms hydrogen bonds with the nicotinamide ribose moiety, thereby lowering the pKa of the TyrOH to promote proton transfer. Hydride transfer is from the Sside of C4 of the nicotinamide to the substrate. AutoDock Vina identified 17 and 15 binding poses for DAUN and DOX, respectively, in the CBR1 binding site and all of them we examined in detail to determine if they could support catalysis. A summary of the docking results is presented in Table 2 where calculated binding affinities and selected distances measuring catalytic efficiency are reported. Among the stable complexes highlighted by AutoDock Vina, arrangements fulfilling the conditions for productive complexes were identified at 8.3 kcal/mol for DAUN and at 7.2 kcal/mol for DOX (pose #3 and #6, respectively, Table 2). These productive poses in the CBR1 active site are displayed in Fig. 3A, and corresponding 2D interaction modes are shown in Fig. 3B,C. In these catalytically productive conformations the C13 carbonyl group is at close distance to the C4 atom of NADPH (3.6 Å and 3.5 Å for DAUN and DOX, respectively), thus allowing for hydride transfer from NADPH to the anthracycline. DOX and DAUN were predicted to bind in the same orientation within the binding site, held in this position by stacking interactions with Trp229 and by additional hydrophobic interactions with the side chain of Met141. This is illustrated in Fig. 3A showing the overlap of DOX and DAUN within CBR1 active site (the RMSD of the heavy atoms is
1.65 Å, whereas the RMSD of the anthracene moiety is 0.17 Å) and how the anthracyclines are sandwiched between residues Trp229 and Met141. As further support for this docking pose, the C13 carbonyl group of both DAUN and DOX is hydrogen bonded to Tyr193 and Ser139 of the catalytic triad (Fig. 3A,B,C). An additional hydrogen bond is established between the anthracycline OH group in C9 and the amide oxygen of the nicotinamide ring of NADPH. On the whole, this binding mode accounts for the reduction of the anthracycline to the corresponding C13-hydroxy metabolite. In the remaining docking poses, the positioning of DAUN and DOX is less favorable for reduction of the C13 carbonyl group. Among them, it is worth mentioning pose #12 for DAUN (AutoDock Vina estimated binding energy ¼ 7.8 kcal/mol) and pose #7 for DOX (AutoDock Vina estimated binding energy ¼ 7.1 kcal/mol). These docked models depict a spatial arrangement between the anthracycline and the cofactor which supports reduction of DAUN and DOX quinone by CBR1. The C4 atom of nicotinamide ring is at 3.9 Å from the quinone C12 atom of the anthracycline, whereas the OH group of Tyr193 is at 4.1 Å and 4.0 Å from the quinone O12 atom of DAUN and DOX, respectively. In these instances also, the anthracene moiety is sandwiched between residues Trp229 and Met141. These results are in good agreement with previous studies by Slupe et al. [13] who manually docked DOX into CBR1 active site and demonstrated that the C13 carbonyl is a preferred substrate for CBR1 in comparison with quinone.
Table 2 Summary of AutoDock Vina results for DAUN and DOX docking to CBR1. Poses in underlined bold represent productive binding modes identified by docking calculations and selected for further analysis. CBR1/DAUN
CBR1/DOX
Pose
DAUN(C13) -NADPH(NC4) distance (Å)
DAUN(O13)TYR(OH) distance (Å)
Predicted affinity (kcal/mol)
Pose
DAUN(C13)NADPH(NC4) distance (Å)
DAUN(O13)TYR(OH) distance (Å)
Predicted affinity (kcal/mol)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
4.3 8.7 3.6 6.7 8.8 3.6 9.7 4.8 8.8 9.5 4.4 7.9 7.6 8.2 12.8 13.5 8.3
2.8 6.8 3.3 5.7 6.9 3.3 9.3 3.1 6.8 11.0 5.2 12.3 12.0 11.6 13.0 14.6 5.8
8.4 8.4 L8.3 8.3 8.2 8.2 8.1 8.0 7.9 7.8 7.8 7.8 7.7 7.6 7.6 7.5 7.4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
8.8 8.7 4.3 6.7 10.0 3.5 8.0 8.6 8.8 7.5 9.3 8.1 5.2 10.8 4.7
6.7 6.6 2.8 4.0 10.3 3.3 12.3 6.3 6.7 11.8 10.8 11.9 4.2 9.9 3.2
7.6 7.6 7.6 7.5 7.4 L7.2 7.1 7.1 7.1 7.1 7.0 6.9 6.8 6.7 6.5
148
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 3. AutoDock Vina-predicted productive poses of DAUN and DOX in CBR1. (A) Shown is the catalytic binding site and the NADPH cofactor in green. The rest of the protein is represented in solid ribbon representation (gold). DOX (red) and DAUN (yellow) are shown as ball and stick models. Predicted hydrogen bonding with Tyr193 and Ser139 are indicated by dashed lines. Orientation was selected for clarity of the p-p interactions predicted between the anthracyclines and Trp229 of CBR1, as well as additional hydrophobic interactions with residue Met141. (B) Schematic representation of the interaction between DAUN and CBR1 produced using the LIGPLOT program [43]. (C) Schematic representation of the interaction between DOX and CBR1 produced using the LIGPLOT program [43]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
In addition, we also examined putative binding modes of epirubicin (EPI) and sabarubicin (MEN 10755) to help elucidate the mechanisms underlying the reduced intramyocardial conversion of EPI and MEN 10755 to secondary alcohol metabolite [21]. Docking simulations of EPI and MEN 10755 into the CBR1 target produced 16 and 15 stable complexes, respectively. In the case of EPI, pose #8 (AutoDock Vina estimated binding energy ¼ 7.1 kcal/mol) overlaps with predicted DOX orientation (RMSD of the heavy atoms ¼ 0.4 Å), whereas pose #4 of MEN 10755 (AutoDock Vina estimated binding energy ¼ 7.7 kcal/mol) is bound in a fashion similar to that seen for DOX, with a RMSD of the anthracene moiety of 5.8 Å. For MEN 10755 binding, its bulky disaccharide moiety compelled the anthracycline to bind in a slightly different orientation relative to that of DOX, with the anthracene moieties coplanar but rotated w15 , such that, in any case, hydride transfer from NADPH is possible (3.6 and 2.9 Å are the catalytically relevant
distances from NADPH and Tyr193, respectively). On the whole docking models cannot suggest a decreased catalytically efficiency of CBR1 toward EPI and MEN 10755, with respect to DOX. 2.2.3. Binding modes of DOX and DAUN with AKR1A1 The 3D structural model of AKR1A1 exhibits a protein fold, an (a/b)8-barrel, characteristic of the AKR superfamily [28]. Reduction of the carbonyl group by AKR1A1 is thought to involve four amino acids: Asp44, Tyr49, Lys79 and His112. Tyr49 is the proton donor, whereas the other amino acids appear to have auxiliary roles in catalysis [28,29]. Asp44/Lys79 are involved in a hydrogen-bonding network while His112 directs the orientation of substrates in the active site pocket. As in the case of CBR1, a minimum energy structure was obtained in which the positional arrangement between the anthracyclines and AKR1A1 can lead to reduction of the C13 carbonyl
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
group. Among the 15 and 13 binding poses identified by AutoDock Vina (Table 3), pose #3 for DAUN (AutoDock Vina estimated binding energy ¼ 9.2 kcal/mol) and pose #2 for DOX (AutoDock Vina estimated binding energy ¼ 9.1 kcal/mol) display the side chain in C9 correctly oriented to undergo the reduction of C13 carbonyl group although the distances from NADPH and Tyr49 are longer than the optimal distances found in CBR1. The binding modes of productive conformations in AKR1A1 are displayed in Fig. 4A and corresponding 2D interaction modes are shown in Fig. 4B,C. The C4 atom of nicotinamide ring is at 4.8 and 4.9 Å from the C13 atom of DAUN and DOX, whereas the OH group of Tyr49 is at 4.3 Å and 4.6 Å from the respectively carbonyl oxygen atoms. Significantly, in both anthracyclines, the C13 oxygen atom is hydrogen bonded to His112 NE2 atom (Fig. 4B,C). Again, DAUN and DOX were predicted to bind in the same orientation within the active site pocket (the RMSD of both heavy atoms and anthracene moiety is 0.24 Å). A number of hydrophobic residues surround the active site, namely Trp21, Ile48, Trp81, Trp113 and Trp219. Both anthracyclines are also held in this position by hydrogen bonding with Ile48, Met301 and Trp21 (Fig. 4A,B,C). This orientation is therefore potentially consistent with the reduction of DAUN and DOX by AKR1A1 but the positioning of the substrate is less favorable for catalysis than found for CBR1. Analogously to CBR1, the location of EPI and MEN 10755 in the binding site of AKR1A1 does not significantly differ from that of DOX. Among the 18 and 17 stable complexes produced by AutoDock Vina for EPI and MEN 10755, respectively, pose #2 of EPI (AutoDock Vina estimated binding energy ¼ 9.1 kcal/mol) overlaps with DOX (RMSD of the heavy atoms ¼ 0.4 Å), whereas pose #13 of MEN 10755 (AutoDock Vina estimated binding energy ¼ 8.6 kcal/mol) adopts a slightly different orientation with respect to DOX (RMSD of the almost coplanar anthracene moieties ¼ 5.9 Å). Catalytically relevant distances from NADPH and Tyr49 (5.0 and 4.9 Å, respectively) do not indicate a decreased catalytically efficiency of AKR1A1 toward EPI and MEN 10755, with respect to DOX. 2.3. Molecular dynamics calculations During the docking runs the ligand was considered as flexible while the protein is fixed in geometry. To account for protein flexibility as well as for the effect of the explicit solvent, the low energy conformations identified by docking calculations, which might undergo the reduction reaction, were submitted to molecular dynamics (MD) calculations. Considering the proteineligand
149
interactions from a dynamic point of view may also allow a more precise evaluation of the arrangements required for catalysis. MD simulations were conducted for all the complexes as well as for the free enzymes. This provided a better picture of the structure and stability of the CBR1- and AKR1A1-anthracycline complexes. Physicochemical parameters such as temperature, potential, and total energy, monitored during the MD simulations, reached stable values after a few hundreds of picoseconds. The Ca-atoms root mean square deviation (RMSD) of all complex and free enzyme simulations is shown in Fig. 5. The RMSD was calculated relative to the initial structure of the production phase during the total 20 ns simulation. From 0 to 500 ps, the RMSD steadily increased as the enzyme structure deviated from its initial structure. After 500 ps, the RMSD was reasonably stable, fluctuating around values of 2.3 0.2 Å with the exception of AKR1A1-DAUN complex (Ca RMSD fluctuating around 2.7 0.1 Å) (Fig. 5). The free enzymes RMSD are very similar to the RMSD of the corresponding enzymesubstrate complexes. We noticed that the overall enzyme structures appeared constant over time in all simulations, except for the segments containing residues 215e226 (loop B) and 298e324 (loop C) of the AKR1A1 enzyme that changed conformation after 100 ps upon DAUN binding. The importance of these findings will be discussed further below. The resulting MD trajectories were also analyzed by monitoring how significant enzymeesubstrate distances varied during the duration of the simulations. The time dependence of the distances between the anthracycline C13 carbon atom and the C4 atom of NADPH, between the anthracycline C13 oxygen atom and the enzyme Tyr193 OH and Ser139 OG atoms, provides insights about the motions of DOX and DAUN at the active site of CBR1 (Fig. 6). With DAUN the distance from NADPH fluctuates around 3.6 0.2 Å all over the simulation time (Fig. 6A), whereas the C13 carbonyl group maintains a strong and stable H-bond with Tyr193 (2.5 0.1 Å) and Ser139 (2.6 0.1 Å) (Fig. 6B and C, respectively). These H-bonds were weakened/broken in the DOX simulation as indicated by the increased separation at 1.5 and 2.2 ns when the distance value of C13 carbonyl group from Ser139 and Tyr193, respectively, has a jump. After the first 2.2 ns, the same distance value between DOX and NADPH increases up to 12 Å, fluctuating around values of 7.3 0.5 Å in the last 5 ns (Fig. 6A). Overall, for the DAUN complex, the distances from the site of catalysis and the number of hydrogen bonds in the docked complex are maintained along the whole 20-ns trajectory. In the case of DOX, the anthracycline during the MD run experienced major adjustments which
Table 3 Summary for AutoDock Vina results for DAUN and DOX docking to AKR1A1. Poses in underlined bold represent productive binding modes identified by docking calculations and selected for further analysis. AKR1A1/DAUN
AKR1A1/DOX
Pose
DAUN(C13)NADPH(NC4) distance ( A)
DAUN(O13)TYR(OH) distance ( A)
Predicted affinity (kcal/mol)
Pose
DAUN(C13)NADPH(NC4) distance ( A)
DAUN(O13)TYR(OH) distance ( A)
Predicted affinity (kcal/mol)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
14.3 10.3 4.8 9.6 14.6 14.7 14.3 10.7 14.8 9.8 6.2 7.6 15.8 9.1 17.9
17.1 11.8 4.3 11.5 14.7 14.3 17.2 12.8 15.0 7.1 7.7 10.1 16.3 6.2 19.3
9.4 9.3 L9.2 8.9 8.8 8.8 8.8 8.8 8.8 8.6 8.3 8.3 7.9 7.9 7.7
1 2 3 4 5 6 7 8 9 10 11 12 13
10.3 4.9 14.3 9.9 9.9 14.8 6.2 14.6 15.8 7.4 7.8 15.7 12.6
12.7 4.6 17.1 11.8 8.9 13.9 7.7 13.4 16.2 9.3 10.5 16.4 11.0
9.4 L9.1 9.0 8.8 8.5 8.3 8.3 8.3 7.9 7.9 7.8 7.8 7.5
150
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 4. AutoDock Vina-predicted productive poses of DAUN and DOX in AKR1A1. (A) Shown is the catalytic binding site and the NADPH cofactor in green. The rest of the protein is represented in solid ribbon representation (gold). DOX (red) and DAUN (yellow) are shown as ball and stick models. Docking of DAUN and DOX is shown to AKR1A1 with the C12 carbonyl oxygen and the 40 -hydroxyl group of the daunosamine moiety hydrogen bonding to Met301 and Ile48, respectively. Hydrophobic residues at the ligand binding site are also shown. (B) Schematic representation of the interaction between DAUN and AKR1A1 produced using the LIGPLOT program [43]. (C) Schematic representation of the interaction between DOX and AKR1A1 produced using the LIGPLOT program [43]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. Ca-RMSDs versus time. RMSD variations of substrate-free CBR1 (green), CBR1/ DAUN (black), CBR1/DOX (red), substrate-free AKR1A1 (brown), AKR1A1/DAUN (blue) and AKR1A1/DOX (yellow) systems with respect to the initial structure as a function of simulation time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
significantly increased the distance of the anthracycline from the CBR1 active site in such a way that its C13 carbonyl group remains close to the catalytic center only for the first 2.2 ns. Analysis of the trajectories reveals that DAUN is more stabilized by additional hydrophobic interactions in the active site, namely by interaction with Met234, whose role in contributing to the preferred substrate specificity for DAUN over that of DOX had be also observed by Slupe et al. [13], and Ala93. The DAUN methyl group in C14 points toward Met234 at the beginning of simulation until 6.5 ns (average distance between the C14 carbon atom and Met234 CE atom ¼ 4.1 0.6), when dihedral angle transitions in the flexible side chain of Met234 are seen. In Fig. 6D the distance between the DAUN C14 atom and the center of mass of Met234 side chain is shown (Fig. 6D). Furthermore, as shown in the inset to Fig. 6D, for more than 70% of the simulation time the distance between C14 atom and the CB atom of Ala93 is less than 4.5 Å, assuming an average value of 3.8 0.2 Å. For comparison, the distance of C14
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 6. Distances accounting for the most relevant anthracycline-CBR1 interactions during the MD simulations. DAUN is black; DOX is gray. (A) Distances between the C13 carbon atom and the C4 atom of NADPH. (B) Distances between the C13 oxygen atom and Tyr193 OH atom. (C) Distances between the C13 oxygen atom and Ser139 OG atom. (D) Distances between the C14 carbon atom and the center of mass of Met234 side chain. Inset, distances between the C14 carbon atom and Ala93 CB atom.
atom from the center of mass of Met234 side chain and Ala93 CB atom in CBR1/DOX simulation is also shown (Fig. 6D). In the case of AKR1A1, the time dependence of the distances between the C13 carbon atom and the C4 atom of NADPH, between the anthracycline C13 oxygen atom and the enzyme Tyr49 OH and His112 NE2 atoms is shown in Fig. 7. The average distance between the NADPH C4 atom and the anthracycline C13 atom was stable throughout the MD simulation of DAUN system with an average value of 4.3 0.4 Å (Fig. 7A). DOX maintains a stable distance from NADPH after 1.2 ns (average distance 4.8 0.4 Å). There are no hydrogen bonds formed between the catalytically critical AKR1A1 Tyr49 and both DAUN and DOX. The distance between the hydroxylic oxygen of Tyr49 and the carbonyl O13 of the anthracyclines takes values around 5.8 0.4 Å and 8.3 1.4 Å in the last 5 ns (DAUN and DOX, respectively, Fig. 7B). It can be observed that at various points during the MD run Tyr49-OH is hydrogen bonded to Lys80. In the case of DAUN the carbonyl O13 atom of the anthracycline receives a hydrogen bond donated by the NE2 atom of His112 for 97% of the time with an average distance of 2.90 0.3 Å (Fig. 7C). By contrast, the results from DOX MD simulation showed a disruption of this hydrogen bonding after 400 ps (Fig. 7C). The distance between DOX carbonyl O13 the NE2 atom of His112 takes values above 7.0 Å between 0.4 and 2.0 ns and after that, it takes values around 3.8 0.7. The decreased distance of DAUN from the AKR1A1 active site may be due to stable hydrophobic interactions between the methyl group in C14 and Trp21. The distance between the DAUN C14 atom and the center of mass of Trp21 indole ring was monitored and shown in comparison with the distance in AKR1A1/
151
Fig. 7. Distances accounting for the most relevant anthracycline-AKR1A1 interactions during the MD simulations. DAUN is black; DOX is gray. (A) Distances between the C13 carbon atom and the C4 atom of NADPH. (B) Distances between the C13 oxygen atom and Tyr49 OH atom. (C) Distances between the C13 oxygen atom and His112 NE2 atom. (D) Distances between the C14 carbon atom and the center of mass of Trp21 residue.
DOX complex (Fig. 7D). After 1.5 ns of simulation, the C14-Trp21 distance becomes significantly smaller in AKR1A1/DAUN system (average distance 4.0 0.4 Å) than in AKR1A1/DOX (average distance 7.4 0.4 Å). Analysis of trajectories indicates that, although fluctuating within its binding site, DAUN is strongly anchored to AKR1A1. As a consequence, flexible loops B and C of AKR1A1, which interact with the anthracene moiety, undergo to significant conformational modifications which are compatible with the higher Ca RMSD observed for AKR1A1/DAUN system (Fig. 5). On the whole, the simulations afford evidence that DAUN is a better substrate for CBR1 and AKR1A1 than DOX and suggest a lower catalytic activity of AKR1A1 toward DAUN and DOX when compared to CBR1. 2.4. Comparison with available experimental data It is important for validation of the modeling to appropriately compare our docking and MD-simulated structures with the available experimental results. Despite the importance of DOX and DAUN reduction, considering that variability in the generation of DAUNol and DOXol may affect both the antitumor effects and the risk of cardiotoxicity, the underlying enzymes are poorly characterized. Reduction of DAUN to DAUNol and of DOX to DOXol is predominantly catalyzed in the cytosol [5]. Kinetic data for reduction of DOX by purified CBR1 and AKR1A1 showed that CBR1 had lower Km and higher turnover number (kcat), and consequently higher catalytic efficiency (kcat/Km), with respect to AKR1A1 ([14], Table 1). This can be explained in terms of the shorter distance
152
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
between the substrate carbonyl group in C13 and the enzyme active site as shown by molecular docking and dynamics simulations and reported in Figs. 3e7. As the hydroxyl group in C14 is replaced by a hydrogen (DAUN) the anthracycline is more tightly bound to the enzymeecoenzyme complex (Figs. 6 and 7) and this is expected to increase the catalytic activity of AKR1A1 which is consistent with the experimental observation that the specific activity of AKR1A1 toward DAUN was higher than that toward DOX ([18], Table 1). It is also interesting to compare the DAUN/CBR1 and DOX/CBR1 binding modeled in the present study with two different experimental measurements which support that DAUN is a better substrate than DOX for CBR1 on the basis of significantly higher kcat/Km values ([13,15], Table 1). Our computational results point out the role of a hydrophobic group in C14 in imparting stability to the productive enzymeesubstrate complex. This stabilization would be impaired for DOX which in the course of the simulation interacts less efficiently with the residues of the catalytic site. In summary, our modeled binding modes of daunorubicin and doxorubicin with cytosolic reductases can reasonably explain all of the previously reported experimental data. 3. Conclusions Anticancer DAUN and DOX reductions are associated with increasing risk of chronic cardiotoxicity. One of the factors underlying the interpatient variation in cardiotoxicity may be linked to cytosolic reductases responsible for metabolizing DOX and DAUN. Several reductase enzymes are involved in the metabolism of these drugs, and distinguishing their individual contributions to the total metabolism is a challenging task. It is with these considerations that we present our molecular modeling findings focusing on the interaction of DOX and DAUN with CBR1 and AKR1A1, the most efficient cytosolic reductases acting on anthracycline antibiotics. We describe a combined study in which homology modeling, docking and MD simulations were employed in order to rationalize the question of how the two different substrates influence the selectivity of CBR1 and AKR1A1. Kinetic studies in vitro indicate that CBR1 and AKR1A1 are less active in metabolizing DOX, with respect to DAUN, by the decrease of kcat/Km, the best kinetic parameter to assess overall substrate efficiency, and, in all cases, CBR1 results a more efficient enzyme (Table 1). The ligandeenzyme interaction models, chosen on the basis of both catalytically productive orientation and docking energy score and further evaluated for dynamic stability of the complexes throughout the MD simulations, are fully consistent with the commonly accepted mechanism of catalysis of CBR1 and AKR1A1 and with previously reported kinetic studies. We also show that the more hydrophobic daunorubicin side chain in C9 may have an important role in the substrate recognition mechanism of these enzymes and that a hydrophobic group in C14 confers stability to the enzymeesubstrate complex. Given that inhibitors of the intramyocardial formation of DOXol and DAUNol metabolites might prove useful to mitigate cardiotoxicity, our results have therefore implications for drug discovery. In future outlooks, the binding models presented here, might serve to design a minimal pharmacophore for the specific enzyme-substrate recognition with the aim to design a novel class of inhibitors of anthracycline reductases.
PP (3-(1-tert-butyl-4-amino-1H-pyrazolo[3,4-d]pyrimidin-3-yl) phenol) molecule. To generate docking targets, inhibitors and crystallographic water molecules were removed from PDB structure. The nicotinamide mononucleotide unit of the cofactor was reduced at C4. The crystallographic structure of AKR1A1 (PDB ID: 2ALR [30]) is missing residues 214e225 and a 3D model structure of AKR1A1 was therefore built by homology modeling. AKR1A1 displays high sequence identity (94%) to porcine aldehyde reductase. Thus, the crystal structure of porcine aldehyde reductase in complex with fidarestat (PDB ID 3H4G [31]) was used as structural template. The sequence alignment produced by ClustalW was used, in conjunction with the program MODELLER [32] implemented in DiscoveryStudio v3.0 (Accelrys Inc.), to generate a set of 50 models of AKR1A1 in complex with fidarestat. The homology model was compared to the 2ALR crystal structure by calculating the rootmean square deviation (RMSD) of the a-carbons in DiscoveryStudio. Only the residues present in the crystal structure were compared to the homology model. RMSD values were also calculated for heavy atoms within 3.5 Å of the cocrystallized fidarestat (residues 21, 48, 49, 112, 113, 298, 299, 300). The stereochemical accuracy of the model was also checked by means of PROCHECK v.3.5.4 [33] and VERIFY3D [34] software. DOX and DAUN structures were extracted from their DNA complexes (PDB ID: 151D [35] and 1D33 [36], respectively) and subjected to 1000 steps of steepest descent energy minimization in vacuo, using GROMACS v. 4.5 [37] and the GROMOS96 43a1 force field [38]. The protonation state 30 -NH3 þ in the daunosamine ring was used. EPI and MEN 10755 were built starting from the structure of DOX using DiscoveryStudio v3.0 (Accelrys Inc.). 4.2. Docking calculations Flexible-ligand docking was carried out using AutoDock Vina 1.1 [39]. This software has been shown to produce, together with an increased efficiency in predicting the experimental binding poses and energies, a 2 orders of magnitude speed-up compared with AutoDock 4 [40]. Vina uses as input files the 3D coordinates of both ligand and receptor, which must be converted into the appropriate format using the AutoDockTools (ADT) program [41]. Polar hydrogens were added by using the Hydrogens module in ADT; after that, Gasteiger partial charges were assigned [42]. For all docking calculations, a grid box size of 30 30 30 points (spacing between the grid points of 1 Å) was used, centered on the mass center of the crystallographic ligand and encompassing all active site atoms. The AutoDock Vina parameter “Exhaustiveness”, which determines how comprehensively the program searches for the lowest energy was set to the value of 500. All other parameters were used as defaults. All rotatable bonds within the ligand were allowed to rotate freely, and the receptor was considered rigid. All computations were carried out on a HP XW8600 workstation running Red Hat Enterprise Linux 5. The resulting docking poses were analyzed with DiscoveryStudio to detect catalytically productive orientation of the substrate within the binding site. Schematic two-dimensional representations of the docking results were produced using LIGPLOT [43]. 4.3. Validation of the docking protocol
4. Experimental section 4.1. Coordinates of proteins and substrates Structure of Protein Data Bank (PDB) entry 1WMA [26] was used for docking simulations of CBR1. The 1WMA structure contains one molecule of CBR1 complexed with an NADPþ and with a hydroxy-
The structures of CBR1 bound to hydroxy-PP (PDB ID: 1WMA) and porcine aldehyde reductase (PDB ID: 3H4G) bound to fidarestat were used for docking validation (correctly identifying the binding site and scoring the receptoreligand interactions). Hydroxy-PP and fidarestat were docked to the targets containing CBR1eNADPH and porcine aldehyde reductase-NADPH, respectively. Ligand docked
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
poses were in close agreement with the crystallographically determined positions. The RMSD of the heavy atoms was only 1.45 Å for hydroxy-PP in CBR1 and 1.12 Å for fidarestat in porcine aldehyde reductase. These control simulations demonstrate that AutoDock Vina is able to generate docking models in the chosen crystal structures and can be used to predict DOX and DAUN binding modes in cytosolic reductases. 4.4. Molecular dynamics calculations Enzymeesubstrate complexes and free enzyme dynamics were conducted using the same method. Simulations and analysis of the trajectories were carried out using GROMACS 4.5 [37] and the GROMOS96 43a1 force field [38]. The topology of the protein and NADPH was created using GROMACS utilities, whereas the topology of the ligand was created using the PRODRG server [44], available at http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg. All starting structures were immersed in a triclinic water box under periodic boundary conditions using a 1.0 nm minimum distance from the protein to the box faces. Each system was then neutralized by Cl and Naþ counterions (CBR1 and AKR1A1 systems, respectively) that were added at random positions to the bulk solvent. The final systems consisted of 2690 and 3258 protein atoms (CBR1 and AKR1A1, respectively) surrounded by w11,000 water molecules. Following steepest descents minimization, the systems were equilibrated under NVT conditions for 100 ps at 300 K, followed by 100 ps under NPT conditions, applying position restraints to the protein and ligand atoms. Finally, all restraints were removed and the molecular dynamics were run for 20 ns at 300 K. The temperature was maintained close to its reference value (300 K) by applying the V-rescale [45] thermostat with a coupling constant of 0.1 ps whereas for the constant pressure (1 atm, isotropic coordinate scaling) the ParrinelloeRahman barostat [46] was used with a relaxation time of 2.0 ps. van der Waals interactions were modeled using 6e12 Lennard-Jones potentials with a 1.4 nm cutoff. Long-range electrostatic interactions were calculated using Particle Mesh Ewald method, with a cutoff for the real space term of 0.9 nm. Covalent bonds were constrained using the LINCS algorithm. The time step employed was 2 fs and the coordinates were saved every 2 ps for analysis which was carried out using the standard GROMACS tools. Acknowledgments Computational resources for MD simulations were provided by CASPUR e Inter-University Consortium for Supercomputing application for University and Research (Standard HPC Grant). References [1] R.B. Weiss, The anthracyclines: will we ever find a better doxorubicin? Semin. Oncol. 19 (1992) 670e686. [2] R.C. Young, R.F. Ozols, C.E. Myers, The anthracycline antineoplastic drugs, N. Engl. J. Med. 305 (1981) 139e153. [3] S.B. Kaye, The multidrug resistance phenotype, Br. J. Cancer 58 (1988) 691e694. [4] L.J. Steinherz, Anthracycline-induced cardiotoxicity, Ann. Intern. Med. 126 (1997) 827. author reply 827e828. [5] G. Minotti, P. Menna, E. Salvatorelli, G. Cairo, L. Gianni, Anthracyclines: molecular advances and pharmacologic developments in antitumor activity and cardiotoxicity, Pharmacol. Rev. 56 (2004) 185e229. [6] R. Zucchi, R. Danesi, Cardiac toxicity of antineoplastic anthracyclines, Curr. Med. Chem. Anticancer Agents 3 (2003) 151e171. [7] R.D. Olson, P.S. Mushlin, D.E. Brenner, S. Fleischer, B.J. Cusack, B.K. Chang, et al., Doxorubicin cardiotoxicity may be caused by its metabolite, doxorubicinol, Proc. Natl. Acad. Sci. U. S. A. 85 (1988) 3585e3589. [8] B.J. Cusack, P.S. Mushlin, L.D. Voulelis, X. Li, R.J. Boucek Jr., R.D. Olson, Daunorubicin-induced cardiac injury in the rabbit: a role for daunorubicinol? Toxicol. Appl. Pharmacol. 118 (1993) 177e185. [9] R.L. Felsted, N.R. Bachur, Mammalian carbonyl reductases, Drug Metab. Rev. 11 (1980) 1e60.
153
[10] A. Mordente, G. Minotti, G.E. Martorana, A. Silvestrini, B. Giardina, E. Meucci, Anthracycline secondary alcohol metabolite formation in human or rabbit heart: biochemical aspects and pharmacologic implications, Biochem. Pharmacol. 66 (2003) 989e998. [11] T. Matsunaga, S. Shintani, A. Hara, Multiplicity of mammalian reductases for xenobiotic carbonyl compounds, Drug Metab. Pharmacokinet. 21 (2006) 1e18. [12] H. Ohara, Y. Miyabe, Y. Deyashiki, K. Matsuura, A. Hara, Reduction of drug ketones by dihydrodiol dehydrogenases, carbonyl reductase and aldehyde reductase of human liver, Biochem. Pharmacol. 50 (1995) 221e227. [13] A. Slupe, B. Williams, C. Larson, L.M. Lee, T. Primbs, A.J. Bruesch, et al., Reduction of 13-deoxydoxorubicin and daunorubicinol anthraquinones by human carbonyl reductase, Cardiovasc. Toxicol. 5 (2005) 365e376. [14] N. Kassner, K. Huse, H.-J. Martin, U. Gödtel-Armbrust, A. Metzger, I. Meineke, et al., Carbonyl reductase 1 is a predominant doxorubicin reductase in the human liver, Drug Metab. Dispos. 36 (2008) 2113e2120. [15] O.S. Bains, M.J. Karkling, T.A. Grigliatti, R.E. Reid, K.W. Riggs, Two nonsynonymous single nucleotide polymorphisms of human carbonyl reductase 1 demonstrate reduced in vitro metabolism of daunorubicin and doxorubicin, Drug Metab. Dispos. 37 (2009) 1107e1114. [16] J.G. Blanco, W.M. Leisenring, V.M. Gonzalez-Covarrubias, T.I. Kawashima, S.M. Davies, M.V. Relling, et al., Genetic polymorphisms in the carbonyl reductase 3 gene CBR3 and the NAD(P)H:quinone oxidoreductase 1 gene NQO1 in patients who developed anthracycline-related congestive heart failure after childhood cancer, Cancer 112 (2008) 2789e2795. [17] R.H. Takahashi, O.S. Bains, T.A. Pfeifer, T.A. Grigliatti, R.E. Reid, K.W. Riggs, Aldo-keto reductase 1C2 fails to metabolize doxorubicin and daunorubicin in vitro, Drug Metab. Dispos. 236 (2008) 991e994. [18] T. O’connor, L.S. Ireland, D.J. Harrison, J.D. Hayes, Major differences exist in the function and tissue-specific expression of human aflatoxin B1 aldehyde reductase and the principal human aldo-keto reductase AKR1 family members, Biochem. J. 343 (1999) 487e504. [19] H.-J. Martin, U. Breyer-Pfaff, V. Wsol, S. Venz, S. Block, E. Maser, Purification and characterization of akr1b10 from human liver: role in carbonyl reduction of xenobiotics, Drug Metab. Dispos. 34 (2006) 464e470. [20] A. Mordente, E. Meucci, A. Silvestrini, G.E. Martorana, B. Giardina, New developments in anthracycline-induced cardiotoxicity, Curr. Med. Chem. 16 (2009) 1656e1672. [21] G. Minotti, S. Licata, A. Saponiero, P. Menna, A.M. Calafiore, G. Di Giammarco, et al., Anthracycline metabolism and toxicity in human myocardium: comparisons between doxorubicin, epirubicin, and a novel disaccharide analogue with a reduced level of formation and [4Fe-4S] reactivity of its secondary alcohol metabolite, Chem. Res. Toxicol. 13 (2000) 1336e1341. [22] W. Fiedler, N. Tchen, J. Bloch, P. Fargeot, R. Sorio, J.B. Vermorken, et al., A study from the EORTC new drug development group: open label phase II study of sabarubicin (MEN-10755) in patients with progressive hormone refractory prostate cancer, Eur. J. Cancer 42 (2006) 200e204. [23] E. Salvatorelli, P. Menna, M. Lusini, E. Covino, G. Minotti, Doxorubicinolone formation and efflux: a salvage pathway against epirubicin accumulation in human heart, J. Pharmacol. Exp. Ther. 329 (2009) 175e184. [24] G. Minotti, S. Recalcati, A. Mordente, G. Liberi, A.M. Calafiore, C. Mancuso, P. Preziosi, G. Cairo, The secondary alcohol metabolite of doxorubicin irreversibly inactivates aconitase/iron regulatory protein-1 in cytosolic fractions from human myocardium, FASEB J. 12 (1998) 541e552. [25] V. Gonzalez-Covarrubias, J.L. Kalabus, J.G. Blanco, Inhibition of polymorphic human carbonyl reductase 1 (CBR1) by the cardioprotectant flavonoid 7monohydroxyethyl rutoside (monoHER), Pharm. Res. 25 (2008) 1730e1734. [26] M. Tanaka, R. Bateman, D. Rauh, E. Vaisberg, S. Ramachandani, C. Zhang, et al., An unbiased cell morphology-based screen for new, biologically active small molecules, PLoS Biol. 3 (2005) e128. [27] U. Oppermann, C. Filling, M. Hult, N. Shafqat, X. Wu, M. Lindh, et al., Shortchain dehydrogenases/reductases (SDR): the 2002 update, Chem. Biol. Interact 143e144 (2003) 247e253. [28] J.M. Jez, M.J. Bennett, B.P. Schlegel, M. Lewis, T.M. Penning, Comparative anatomy of the aldo-keto reductase superfamily, Biochem. J. 326 (Pt 3) (1997) 625e636. [29] B.P. Schlegel, J.M. Jez, T.M. Penning, Mutagenesis of 3 alpha-hydroxysteroid dehydrogenase reveals a «push-pull» mechanism for proton transfer in aldo-keto reductases, Biochemistry 37 (1998) 3538e3548. [30] O. El-Kabbani, N.C. Green, G. Lin, M. Carson, S.V. Narayana, K.M. Moore, et al., Structures of human and porcine aldehyde reductase: an enzyme implicated in diabetic complications, Acta Crystallogr. D Biol. Crystallogr. 50 (1994) 859e868. [31] O. El-Kabbani, V. Carbone, C. Darmanin, M. Oka, A. Mitschler, A. Podjarny, et al., Structure of aldehyde reductase holoenzyme in complex with the potent aldose reductase inhibitor fidarestat: implications for inhibitor binding and selectivity, J. Med. Chem. 48 (2005) 5536e5542. [32] A. Sali, T.L. Blundell, Comparative protein modelling by satisfaction of spatial restraints, J. Mol. Biol. 234 (1993) 779e815. [33] A.L. Morris, M.W. MacArthur, E.G. Hutchinson, J.M. Thornton, Stereochemical quality of protein structure coordinates, Proteins 12 (1992) 345e364. [34] R. Lüthy, J.U. Bowie, D. Eisenberg, Assessment of protein models with threedimensional profiles, Nature 356 (1992) 83e85. [35] L.A. Lipscomb, M.E. Peek, F.X. Zhou, J.A. Bertrand, D. VanDerveer, L.D. Williams, Water ring structure at DNA interfaces: hydration and dynamics of DNAanthracycline complexes, Biochemistry 33 (1994) 3649e3659.
154
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
[36] A.H. Wang, Y.G. Gao, Y.C. Liaw, Y.K. Li, Formaldehyde cross-links daunorubicin and DNA efficiently: HPLC and X-ray diffraction studies, Biochemistry 30 (1991) 3812e3815. [37] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theory Comput. 4 (2008) 435e447. [38] W.R.P. Scott, P.H. Hünenberger, I.G. Tironi, A.E. Mark, S.R. Billeter, J. Fennen, et al., The GROMOS biomolecular simulation program package, J. Phys. Chem. A 103 (1999) 3596e3607. [39] O. Trott, A.J. Olson, AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading, J. Comput. Chem. 31 (2010) 455e461. [40] R. Huey, G.M. Morris, A.J. Olson, D.S. Goodsell, A semiempirical free energy force field with charge-based desolvation, J. Comput. Chem. 28 (2007) 1145e1152.
[41] M.F. Sanner, Python: a programming language for software integration and development, J. Mol. Graph. Model. 17 (1999) 57e61. [42] J. Gasteiger, Iterative partial equalization of orbital electronegativity-a rapid access to atomic charges, Tetrahedron 36 (1980) 3219e3228. [43] A.C. Wallace, R.A. Laskowski, J.M. Thornton, LIGPLOT e a program to generate schematic diagrams of protein ligand interactions, Protein Eng. 8 (1995) 127e134. [44] A.W. Schüttelkopf, D.M.F. van Aalten, PRODRG: a tool for high-throughput crystallography of protein-ligand complexes, Acta Crystallogr. D Biol. Crystallogr. 60 (2204) 1355e1363. [45] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (2007) 014101. [46] M. Parrinello, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182.
Contents lists available at SciVerse ScienceDirect
European Journal of Medicinal Chemistry journal homepage: http://www.elsevier.com/locate/ejmech
Original article
Understanding the binding of daunorubicin and doxorubicin to NADPH-dependent cytosolic reductases by computational methods Davide Pirolli a, Bruno Giardina a, Alvaro Mordente a, Silvana Ficarra b, Maria Cristina De Rosa c, * a
Istituto di Biochimica e Biochimica Clinica, Università Cattolica del Sacro Cuore, Largo F. Vito 1, 00168 Rome, Italy Dipartimento di Chimica Organica e Biologica, Università di Messina, V. le Ferdinando Stagno d’Alcontres 31, 98166 Messina, Italy c CNR-Istituto di Chimica del Riconoscimento Molecolare, c/o Istituto di Biochimica e Biochimica Clinica, Università Cattolica del Sacro Cuore, Largo F. Vito 1, I-00168 Rome, Italy b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 May 2012 Received in revised form 2 August 2012 Accepted 16 August 2012 Available online 24 August 2012
The anthracycline anticancer agents daunorubicin (DAUN) and doxorubicin (DOX) are reduced by different NADPH-dependent cytosolic reductases into their corresponding alcohol metabolites daunorubicinol (DAUNol) and doxorubicinol (DOXol), which have been implicated in the development of chronic cardiomyopathy. To better understand the individual importance of each enzyme in the reduction and to provide deeper insight into the binding at atomic level we performed molecular docking and dynamics simulations of DAUN and DOX into the active sites of human carbonyl reductase 1 (CBR1) and human aldehyde reductase (AKR1A1). Such simulations evidenced a different behavior between the reductases with respect to DAUN and DOX suggesting major contribution of CBR1 in the reduction. The results are in agreement with available experimental data and for each enzyme and anthracycline pair provided the identification of key residues involved in the interactions. The structural models that we have derived could serve as a useful tool for structure-guided drug design studies. Ó 2012 Elsevier Masson SAS. All rights reserved.
Keywords: Anthracycline Cardiotoxicity Molecular docking Molecular dynamics
1. Introduction The anthracycline antibiotics doxorubicin (DOX) and its congener, daunorubicin (DAUN) are among the most effective chemotherapeutic agents in the treatment of numerous malignancies [1]. Since their first clinical trials these drugs have been extensively studied alone and in combination, and they now have a major role in the effective treatment of acute leukemia, nonHodgkin lymphomas, breast cancer, Hodgkin’s disease, and sarcomas [2]. However, two major problems limit their clinical application: drug resistance [3] and cardiotoxicity [4]. This last complication represents a major obstacle to prolonged use of the drug and/or cumulative dose exceeding 550 mg/m2 body surface area [5]. A number of mechanisms have been proposed to account for the anthracycline-induced cardiotoxicity [6]. Among them, a prevailing hypothesis implicates the C-13 hydroxy metabolite, doxorubicinol (DOXol) and daunorubicinol (DAUNol), compounds that lack the antiproliferative activity of DOX and DAUN, as a major component in the development of chronic cardiomyopathy [7]. The
* Corresponding author. E-mail addresses: [email protected] (D. Pirolli), bgiardina@ rm.unicatt.it (B. Giardina), [email protected] (A. Mordente), sficarra@ isengard.unime.it (S. Ficarra), [email protected] (M.C. De Rosa). 0223-5234/$ e see front matter Ó 2012 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ejmech.2012.08.023
accumulation of the metabolite in the heart has been described suggesting that the C-13 hydroxy metabolite might be transported and selectively trapped by cardiac cells [8]. Several cytosolic reductases of the short chain dehydrogenase class (SDR) and of the aldoeketo reductase class (AKR) metabolize DOX to DOXol and DAUN to DAUNol (Fig. 1) implicating these enzymes in the development of cardiotoxicity [9,10]. Enzymes that have actually been described to reduce anthracyclines to C13-hydroxymetabolites are carbonyl reductase 1 (CBR1) [11e15], carbonyl reductase 3 (CBR3) [16], aldo-keto reductase (AKR) 1A1 [12,14,17], AKR1B1 [14,18,19], AKR1B10 [14,18,19], AKR1C3 [14] and AKR1C4 [14]. Divergent results are reported on the ability of AKR1C2 to metabolize doxorubicin and daunorubicin in vitro [12,17]. At present, CBR1 and AKR1A1 are considered the major anthracycline-metabolizing enzymes, at least in human heart [20] and indications about DOX and DAUN reduction obtained from experimental studies are reported in Table 1. Much effort has been devoted toward developing novel anthracyclines showing increased cardiac tolerability with little success. A decrease in the level of alcohol metabolite formation was observed with epirubicin (EPI), the C-40 epimer of DOX, and sabarubicin (MEN 10755), a DOX analog lacking the methoxy group at C-4 and characterized with a disaccharide in the place of daunosamine in C-7 [21]. However these anthracyclines are still characterized by a low but significant cardiotoxicity that makes necessary further investigations toward a less toxic anticancer
146
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 1. Enzymatic reduction of anthracyclines to C-13 secondary alcohol metabolites. Schematic representation of two-electron reduction of the side chain carbonyl group of DOX and DAUN resulting in DOXol and DAUNol, respectively.
anthracycline [22,23]. The evaluation of efficient inhibitors of anthracycline reductases has been also suggested as one potential strategy to optimize cardioprotection [24,25], but structure-based design of specific inhibitors is lacking. On the whole, despite the large amount of available biochemical data, there is no structural information at the atomic level regarding these interactions, how cytosolic reductases bind the anthracyclines within the protein cavity and which residues are involved in the binding process. In this study we have used computational methods, based on molecular docking studies and molecular dynamics (MD) simulation, to explore these features in some detail and to provide useful insights into the future rational design of specific inhibitors of anthracycline reductases activity as well as of novel, less cardiotoxic anthracyclines. Our computational studies correlated with previous knowledge in the field and provide a plausible binding model of DAUN and DOX within the active site of CBR1 and AKR1A1. 2. Results and discussion
2.2. Docking calculations
In carrying out the analysis of DOX and DAUN binding in cytosolic reductases, substrates were first docked into the enzyme active site. In this step the substrate was flexible while the enzyme was fixed in geometry. To account for protein flexibility and to allow the geometry to relax, in a second step the low energy docked conformations were submitted to molecular dynamics simulations. Starting coordinates for CBR1 were obtained from the Protein Data Bank while for AKR1A1 atomic coordinates were generated by homology modeling. 2.1. Construction of AKR1A1 model The structure of human aldehyde reductase at the resolution of 2.48 Å (PDB ID: 2ALR) does not completely resolve the polypeptide chain. In fact, residues 214e225 belonging to the extremely flexible Table 1 Experimental measurements, obtained from the literature, for the reaction of CBR1 and AKR1A1 with the substrates DOX and DAUN. Kinetic parameter
kcat(s1) Km(mM) Vmax/Km (U/mg/M) kcat(min1) Km(mM) kcat(s1) Km(mM) Specific activitya a
CBR1
AKR1A1
Ref
DAUN
DOX
DAUN
DOX
1.27 0.03 108 7 2140 e e 1.94 51 13 e
0.27 0.02 90 16 e 0.62 167 0.21 287 57 e
e e 312 e e e e 120 3
e e e 0.04 247 e e 11 1
nmol/min per mg of protein.
loop B of aldehyde reductase are missing in the 2ALR structure. A 3D structure model of AKR1A was then obtained using homology method. The appropriate template was chosen based on sequence similarity, residue completeness, and crystal resolution. The target and the template, porcine aldehyde reductase (PDB ID: 3H4G), share 94% of their sequences which guarantees the accuracy of the prediction. Among the conformations generated by MODELLER one was selected on the basis of the lowest MODELLER objective function. The Ca root-mean square deviation (RMSD) of the homology model to the template structure was 0.24 Å. The homology model was also compared to the human aldehyde reductase crystal structure (2ALR) comparing only the residues present in the crystal structure. RMSD values were also calculated for heavy atoms within 3.5 Å of the co-crystallized fidarestat (residues 21, 48, 49, 112, 113, 298, 299, 300). The Ca RMSD for the model and the AKR1A crystal structure (2ALR) was 0.9 Å whereas the RMSD of fidarestat binding site residues, within 3.5 Å of the ligand was 1.2 Å. The stereochemical quality of the selected model was examined using the program PROCHECK. The distribution of the F and J angles in the Ramachandran plot is well within the allowed regions. The majority (94%) of the residues are in the most favored regions of the plot; 6% are in the additional allowed regions; no residues lie in disallowed regions. The final structure was further checked by Verify3D, as implemented in the Verify3D Structure Evaluation Server (http://nihserver.mbi.ucla.edu/Verify_ 3D/), a program that analyzes the compatibility of an atomic model (3D) with its own amino acid sequence (1D). As shown in Fig. 2, the VERIFY-3D scores of our model are always positive and similar to those of the template protein and those of human aldehyde reductase. Taken together, these data suggest that built 3D model of AKR1A1 was of reasonable quality compared to 3H4G and 2ALR crystal structures.
[13] [13] [12] [14] [14] [15] [15] [18]
Cytosolic reductases were shown to follow a sequential ordered mechanism in which NADPH binds to the enzyme first and NADPþ is released last after the alcohol product. The reaction mechanism involves a hydride transfer from NADPH to the carbonyl carbon of the substrate and a proton transfer from a Tyr residue at the bottom of the active site cavity, which serves as the general acid catalyst, to the forming alcohol group. Additional active site residues form a proton relay system and stabilize the orientation of the substrate. The conformations from the docking experiments were therefore analyzed for a catalytically productive orientation of the substrate within the binding site. Two criteria were used to identify productive conformations: 1) the distance between the C13 anthracycline carbon atom and the C4 atom of NADPH should be short enough to allow efficient hydride transfer; 2) the oxygen of C13 carbonyl should be located in proximity to the Tyr residue which acts as proton donor during the catalysis. 2.2.1. Test of the method The structures of CBR1 bound to hydroxy-PP and porcine aldehyde reductase bound to fidarestat were chosen to test the docking protocol. The results of docking validation indicate that the predicted binding conformations determined by AutoDock Vina match well with those of the co-crystallized ligand. 2.2.2. Binding modes of DAUN and DOX with CBR1 The overall structure of CBR1 retains the characteristic a/ b folding pattern with a Rossmann-fold of the SDR family and three residues (Ser139, Tyr193 and Lys197) were proposed as the catalytic triad in this short-chain dehydrogenase/reductase enzyme [26,27]. Accordingly, Tyr193 functions as the catalytic acid/base,
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
147
Fig. 2. Average Verify3D score profile. Plots obtained by the Verify-3D program, with a sliding 21-residue window, using the coordinates of AKR1A1 homology model (black), the template protein 3H4G (black dashed) and the crystallographic aldehyde reductase 2ALR (gray). Note that residues 214e224 were unresolved in the 2ALR structure and are therefore missing in the corresponding plot.
Ser139 stabilizes the substrate by forming interactions to the substrate carbonyl, and Lys197 forms hydrogen bonds with the nicotinamide ribose moiety, thereby lowering the pKa of the TyrOH to promote proton transfer. Hydride transfer is from the Sside of C4 of the nicotinamide to the substrate. AutoDock Vina identified 17 and 15 binding poses for DAUN and DOX, respectively, in the CBR1 binding site and all of them we examined in detail to determine if they could support catalysis. A summary of the docking results is presented in Table 2 where calculated binding affinities and selected distances measuring catalytic efficiency are reported. Among the stable complexes highlighted by AutoDock Vina, arrangements fulfilling the conditions for productive complexes were identified at 8.3 kcal/mol for DAUN and at 7.2 kcal/mol for DOX (pose #3 and #6, respectively, Table 2). These productive poses in the CBR1 active site are displayed in Fig. 3A, and corresponding 2D interaction modes are shown in Fig. 3B,C. In these catalytically productive conformations the C13 carbonyl group is at close distance to the C4 atom of NADPH (3.6 Å and 3.5 Å for DAUN and DOX, respectively), thus allowing for hydride transfer from NADPH to the anthracycline. DOX and DAUN were predicted to bind in the same orientation within the binding site, held in this position by stacking interactions with Trp229 and by additional hydrophobic interactions with the side chain of Met141. This is illustrated in Fig. 3A showing the overlap of DOX and DAUN within CBR1 active site (the RMSD of the heavy atoms is
1.65 Å, whereas the RMSD of the anthracene moiety is 0.17 Å) and how the anthracyclines are sandwiched between residues Trp229 and Met141. As further support for this docking pose, the C13 carbonyl group of both DAUN and DOX is hydrogen bonded to Tyr193 and Ser139 of the catalytic triad (Fig. 3A,B,C). An additional hydrogen bond is established between the anthracycline OH group in C9 and the amide oxygen of the nicotinamide ring of NADPH. On the whole, this binding mode accounts for the reduction of the anthracycline to the corresponding C13-hydroxy metabolite. In the remaining docking poses, the positioning of DAUN and DOX is less favorable for reduction of the C13 carbonyl group. Among them, it is worth mentioning pose #12 for DAUN (AutoDock Vina estimated binding energy ¼ 7.8 kcal/mol) and pose #7 for DOX (AutoDock Vina estimated binding energy ¼ 7.1 kcal/mol). These docked models depict a spatial arrangement between the anthracycline and the cofactor which supports reduction of DAUN and DOX quinone by CBR1. The C4 atom of nicotinamide ring is at 3.9 Å from the quinone C12 atom of the anthracycline, whereas the OH group of Tyr193 is at 4.1 Å and 4.0 Å from the quinone O12 atom of DAUN and DOX, respectively. In these instances also, the anthracene moiety is sandwiched between residues Trp229 and Met141. These results are in good agreement with previous studies by Slupe et al. [13] who manually docked DOX into CBR1 active site and demonstrated that the C13 carbonyl is a preferred substrate for CBR1 in comparison with quinone.
Table 2 Summary of AutoDock Vina results for DAUN and DOX docking to CBR1. Poses in underlined bold represent productive binding modes identified by docking calculations and selected for further analysis. CBR1/DAUN
CBR1/DOX
Pose
DAUN(C13) -NADPH(NC4) distance (Å)
DAUN(O13)TYR(OH) distance (Å)
Predicted affinity (kcal/mol)
Pose
DAUN(C13)NADPH(NC4) distance (Å)
DAUN(O13)TYR(OH) distance (Å)
Predicted affinity (kcal/mol)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
4.3 8.7 3.6 6.7 8.8 3.6 9.7 4.8 8.8 9.5 4.4 7.9 7.6 8.2 12.8 13.5 8.3
2.8 6.8 3.3 5.7 6.9 3.3 9.3 3.1 6.8 11.0 5.2 12.3 12.0 11.6 13.0 14.6 5.8
8.4 8.4 L8.3 8.3 8.2 8.2 8.1 8.0 7.9 7.8 7.8 7.8 7.7 7.6 7.6 7.5 7.4
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
8.8 8.7 4.3 6.7 10.0 3.5 8.0 8.6 8.8 7.5 9.3 8.1 5.2 10.8 4.7
6.7 6.6 2.8 4.0 10.3 3.3 12.3 6.3 6.7 11.8 10.8 11.9 4.2 9.9 3.2
7.6 7.6 7.6 7.5 7.4 L7.2 7.1 7.1 7.1 7.1 7.0 6.9 6.8 6.7 6.5
148
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 3. AutoDock Vina-predicted productive poses of DAUN and DOX in CBR1. (A) Shown is the catalytic binding site and the NADPH cofactor in green. The rest of the protein is represented in solid ribbon representation (gold). DOX (red) and DAUN (yellow) are shown as ball and stick models. Predicted hydrogen bonding with Tyr193 and Ser139 are indicated by dashed lines. Orientation was selected for clarity of the p-p interactions predicted between the anthracyclines and Trp229 of CBR1, as well as additional hydrophobic interactions with residue Met141. (B) Schematic representation of the interaction between DAUN and CBR1 produced using the LIGPLOT program [43]. (C) Schematic representation of the interaction between DOX and CBR1 produced using the LIGPLOT program [43]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
In addition, we also examined putative binding modes of epirubicin (EPI) and sabarubicin (MEN 10755) to help elucidate the mechanisms underlying the reduced intramyocardial conversion of EPI and MEN 10755 to secondary alcohol metabolite [21]. Docking simulations of EPI and MEN 10755 into the CBR1 target produced 16 and 15 stable complexes, respectively. In the case of EPI, pose #8 (AutoDock Vina estimated binding energy ¼ 7.1 kcal/mol) overlaps with predicted DOX orientation (RMSD of the heavy atoms ¼ 0.4 Å), whereas pose #4 of MEN 10755 (AutoDock Vina estimated binding energy ¼ 7.7 kcal/mol) is bound in a fashion similar to that seen for DOX, with a RMSD of the anthracene moiety of 5.8 Å. For MEN 10755 binding, its bulky disaccharide moiety compelled the anthracycline to bind in a slightly different orientation relative to that of DOX, with the anthracene moieties coplanar but rotated w15 , such that, in any case, hydride transfer from NADPH is possible (3.6 and 2.9 Å are the catalytically relevant
distances from NADPH and Tyr193, respectively). On the whole docking models cannot suggest a decreased catalytically efficiency of CBR1 toward EPI and MEN 10755, with respect to DOX. 2.2.3. Binding modes of DOX and DAUN with AKR1A1 The 3D structural model of AKR1A1 exhibits a protein fold, an (a/b)8-barrel, characteristic of the AKR superfamily [28]. Reduction of the carbonyl group by AKR1A1 is thought to involve four amino acids: Asp44, Tyr49, Lys79 and His112. Tyr49 is the proton donor, whereas the other amino acids appear to have auxiliary roles in catalysis [28,29]. Asp44/Lys79 are involved in a hydrogen-bonding network while His112 directs the orientation of substrates in the active site pocket. As in the case of CBR1, a minimum energy structure was obtained in which the positional arrangement between the anthracyclines and AKR1A1 can lead to reduction of the C13 carbonyl
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
group. Among the 15 and 13 binding poses identified by AutoDock Vina (Table 3), pose #3 for DAUN (AutoDock Vina estimated binding energy ¼ 9.2 kcal/mol) and pose #2 for DOX (AutoDock Vina estimated binding energy ¼ 9.1 kcal/mol) display the side chain in C9 correctly oriented to undergo the reduction of C13 carbonyl group although the distances from NADPH and Tyr49 are longer than the optimal distances found in CBR1. The binding modes of productive conformations in AKR1A1 are displayed in Fig. 4A and corresponding 2D interaction modes are shown in Fig. 4B,C. The C4 atom of nicotinamide ring is at 4.8 and 4.9 Å from the C13 atom of DAUN and DOX, whereas the OH group of Tyr49 is at 4.3 Å and 4.6 Å from the respectively carbonyl oxygen atoms. Significantly, in both anthracyclines, the C13 oxygen atom is hydrogen bonded to His112 NE2 atom (Fig. 4B,C). Again, DAUN and DOX were predicted to bind in the same orientation within the active site pocket (the RMSD of both heavy atoms and anthracene moiety is 0.24 Å). A number of hydrophobic residues surround the active site, namely Trp21, Ile48, Trp81, Trp113 and Trp219. Both anthracyclines are also held in this position by hydrogen bonding with Ile48, Met301 and Trp21 (Fig. 4A,B,C). This orientation is therefore potentially consistent with the reduction of DAUN and DOX by AKR1A1 but the positioning of the substrate is less favorable for catalysis than found for CBR1. Analogously to CBR1, the location of EPI and MEN 10755 in the binding site of AKR1A1 does not significantly differ from that of DOX. Among the 18 and 17 stable complexes produced by AutoDock Vina for EPI and MEN 10755, respectively, pose #2 of EPI (AutoDock Vina estimated binding energy ¼ 9.1 kcal/mol) overlaps with DOX (RMSD of the heavy atoms ¼ 0.4 Å), whereas pose #13 of MEN 10755 (AutoDock Vina estimated binding energy ¼ 8.6 kcal/mol) adopts a slightly different orientation with respect to DOX (RMSD of the almost coplanar anthracene moieties ¼ 5.9 Å). Catalytically relevant distances from NADPH and Tyr49 (5.0 and 4.9 Å, respectively) do not indicate a decreased catalytically efficiency of AKR1A1 toward EPI and MEN 10755, with respect to DOX. 2.3. Molecular dynamics calculations During the docking runs the ligand was considered as flexible while the protein is fixed in geometry. To account for protein flexibility as well as for the effect of the explicit solvent, the low energy conformations identified by docking calculations, which might undergo the reduction reaction, were submitted to molecular dynamics (MD) calculations. Considering the proteineligand
149
interactions from a dynamic point of view may also allow a more precise evaluation of the arrangements required for catalysis. MD simulations were conducted for all the complexes as well as for the free enzymes. This provided a better picture of the structure and stability of the CBR1- and AKR1A1-anthracycline complexes. Physicochemical parameters such as temperature, potential, and total energy, monitored during the MD simulations, reached stable values after a few hundreds of picoseconds. The Ca-atoms root mean square deviation (RMSD) of all complex and free enzyme simulations is shown in Fig. 5. The RMSD was calculated relative to the initial structure of the production phase during the total 20 ns simulation. From 0 to 500 ps, the RMSD steadily increased as the enzyme structure deviated from its initial structure. After 500 ps, the RMSD was reasonably stable, fluctuating around values of 2.3 0.2 Å with the exception of AKR1A1-DAUN complex (Ca RMSD fluctuating around 2.7 0.1 Å) (Fig. 5). The free enzymes RMSD are very similar to the RMSD of the corresponding enzymesubstrate complexes. We noticed that the overall enzyme structures appeared constant over time in all simulations, except for the segments containing residues 215e226 (loop B) and 298e324 (loop C) of the AKR1A1 enzyme that changed conformation after 100 ps upon DAUN binding. The importance of these findings will be discussed further below. The resulting MD trajectories were also analyzed by monitoring how significant enzymeesubstrate distances varied during the duration of the simulations. The time dependence of the distances between the anthracycline C13 carbon atom and the C4 atom of NADPH, between the anthracycline C13 oxygen atom and the enzyme Tyr193 OH and Ser139 OG atoms, provides insights about the motions of DOX and DAUN at the active site of CBR1 (Fig. 6). With DAUN the distance from NADPH fluctuates around 3.6 0.2 Å all over the simulation time (Fig. 6A), whereas the C13 carbonyl group maintains a strong and stable H-bond with Tyr193 (2.5 0.1 Å) and Ser139 (2.6 0.1 Å) (Fig. 6B and C, respectively). These H-bonds were weakened/broken in the DOX simulation as indicated by the increased separation at 1.5 and 2.2 ns when the distance value of C13 carbonyl group from Ser139 and Tyr193, respectively, has a jump. After the first 2.2 ns, the same distance value between DOX and NADPH increases up to 12 Å, fluctuating around values of 7.3 0.5 Å in the last 5 ns (Fig. 6A). Overall, for the DAUN complex, the distances from the site of catalysis and the number of hydrogen bonds in the docked complex are maintained along the whole 20-ns trajectory. In the case of DOX, the anthracycline during the MD run experienced major adjustments which
Table 3 Summary for AutoDock Vina results for DAUN and DOX docking to AKR1A1. Poses in underlined bold represent productive binding modes identified by docking calculations and selected for further analysis. AKR1A1/DAUN
AKR1A1/DOX
Pose
DAUN(C13)NADPH(NC4) distance ( A)
DAUN(O13)TYR(OH) distance ( A)
Predicted affinity (kcal/mol)
Pose
DAUN(C13)NADPH(NC4) distance ( A)
DAUN(O13)TYR(OH) distance ( A)
Predicted affinity (kcal/mol)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
14.3 10.3 4.8 9.6 14.6 14.7 14.3 10.7 14.8 9.8 6.2 7.6 15.8 9.1 17.9
17.1 11.8 4.3 11.5 14.7 14.3 17.2 12.8 15.0 7.1 7.7 10.1 16.3 6.2 19.3
9.4 9.3 L9.2 8.9 8.8 8.8 8.8 8.8 8.8 8.6 8.3 8.3 7.9 7.9 7.7
1 2 3 4 5 6 7 8 9 10 11 12 13
10.3 4.9 14.3 9.9 9.9 14.8 6.2 14.6 15.8 7.4 7.8 15.7 12.6
12.7 4.6 17.1 11.8 8.9 13.9 7.7 13.4 16.2 9.3 10.5 16.4 11.0
9.4 L9.1 9.0 8.8 8.5 8.3 8.3 8.3 7.9 7.9 7.8 7.8 7.5
150
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 4. AutoDock Vina-predicted productive poses of DAUN and DOX in AKR1A1. (A) Shown is the catalytic binding site and the NADPH cofactor in green. The rest of the protein is represented in solid ribbon representation (gold). DOX (red) and DAUN (yellow) are shown as ball and stick models. Docking of DAUN and DOX is shown to AKR1A1 with the C12 carbonyl oxygen and the 40 -hydroxyl group of the daunosamine moiety hydrogen bonding to Met301 and Ile48, respectively. Hydrophobic residues at the ligand binding site are also shown. (B) Schematic representation of the interaction between DAUN and AKR1A1 produced using the LIGPLOT program [43]. (C) Schematic representation of the interaction between DOX and AKR1A1 produced using the LIGPLOT program [43]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. Ca-RMSDs versus time. RMSD variations of substrate-free CBR1 (green), CBR1/ DAUN (black), CBR1/DOX (red), substrate-free AKR1A1 (brown), AKR1A1/DAUN (blue) and AKR1A1/DOX (yellow) systems with respect to the initial structure as a function of simulation time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
significantly increased the distance of the anthracycline from the CBR1 active site in such a way that its C13 carbonyl group remains close to the catalytic center only for the first 2.2 ns. Analysis of the trajectories reveals that DAUN is more stabilized by additional hydrophobic interactions in the active site, namely by interaction with Met234, whose role in contributing to the preferred substrate specificity for DAUN over that of DOX had be also observed by Slupe et al. [13], and Ala93. The DAUN methyl group in C14 points toward Met234 at the beginning of simulation until 6.5 ns (average distance between the C14 carbon atom and Met234 CE atom ¼ 4.1 0.6), when dihedral angle transitions in the flexible side chain of Met234 are seen. In Fig. 6D the distance between the DAUN C14 atom and the center of mass of Met234 side chain is shown (Fig. 6D). Furthermore, as shown in the inset to Fig. 6D, for more than 70% of the simulation time the distance between C14 atom and the CB atom of Ala93 is less than 4.5 Å, assuming an average value of 3.8 0.2 Å. For comparison, the distance of C14
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
Fig. 6. Distances accounting for the most relevant anthracycline-CBR1 interactions during the MD simulations. DAUN is black; DOX is gray. (A) Distances between the C13 carbon atom and the C4 atom of NADPH. (B) Distances between the C13 oxygen atom and Tyr193 OH atom. (C) Distances between the C13 oxygen atom and Ser139 OG atom. (D) Distances between the C14 carbon atom and the center of mass of Met234 side chain. Inset, distances between the C14 carbon atom and Ala93 CB atom.
atom from the center of mass of Met234 side chain and Ala93 CB atom in CBR1/DOX simulation is also shown (Fig. 6D). In the case of AKR1A1, the time dependence of the distances between the C13 carbon atom and the C4 atom of NADPH, between the anthracycline C13 oxygen atom and the enzyme Tyr49 OH and His112 NE2 atoms is shown in Fig. 7. The average distance between the NADPH C4 atom and the anthracycline C13 atom was stable throughout the MD simulation of DAUN system with an average value of 4.3 0.4 Å (Fig. 7A). DOX maintains a stable distance from NADPH after 1.2 ns (average distance 4.8 0.4 Å). There are no hydrogen bonds formed between the catalytically critical AKR1A1 Tyr49 and both DAUN and DOX. The distance between the hydroxylic oxygen of Tyr49 and the carbonyl O13 of the anthracyclines takes values around 5.8 0.4 Å and 8.3 1.4 Å in the last 5 ns (DAUN and DOX, respectively, Fig. 7B). It can be observed that at various points during the MD run Tyr49-OH is hydrogen bonded to Lys80. In the case of DAUN the carbonyl O13 atom of the anthracycline receives a hydrogen bond donated by the NE2 atom of His112 for 97% of the time with an average distance of 2.90 0.3 Å (Fig. 7C). By contrast, the results from DOX MD simulation showed a disruption of this hydrogen bonding after 400 ps (Fig. 7C). The distance between DOX carbonyl O13 the NE2 atom of His112 takes values above 7.0 Å between 0.4 and 2.0 ns and after that, it takes values around 3.8 0.7. The decreased distance of DAUN from the AKR1A1 active site may be due to stable hydrophobic interactions between the methyl group in C14 and Trp21. The distance between the DAUN C14 atom and the center of mass of Trp21 indole ring was monitored and shown in comparison with the distance in AKR1A1/
151
Fig. 7. Distances accounting for the most relevant anthracycline-AKR1A1 interactions during the MD simulations. DAUN is black; DOX is gray. (A) Distances between the C13 carbon atom and the C4 atom of NADPH. (B) Distances between the C13 oxygen atom and Tyr49 OH atom. (C) Distances between the C13 oxygen atom and His112 NE2 atom. (D) Distances between the C14 carbon atom and the center of mass of Trp21 residue.
DOX complex (Fig. 7D). After 1.5 ns of simulation, the C14-Trp21 distance becomes significantly smaller in AKR1A1/DAUN system (average distance 4.0 0.4 Å) than in AKR1A1/DOX (average distance 7.4 0.4 Å). Analysis of trajectories indicates that, although fluctuating within its binding site, DAUN is strongly anchored to AKR1A1. As a consequence, flexible loops B and C of AKR1A1, which interact with the anthracene moiety, undergo to significant conformational modifications which are compatible with the higher Ca RMSD observed for AKR1A1/DAUN system (Fig. 5). On the whole, the simulations afford evidence that DAUN is a better substrate for CBR1 and AKR1A1 than DOX and suggest a lower catalytic activity of AKR1A1 toward DAUN and DOX when compared to CBR1. 2.4. Comparison with available experimental data It is important for validation of the modeling to appropriately compare our docking and MD-simulated structures with the available experimental results. Despite the importance of DOX and DAUN reduction, considering that variability in the generation of DAUNol and DOXol may affect both the antitumor effects and the risk of cardiotoxicity, the underlying enzymes are poorly characterized. Reduction of DAUN to DAUNol and of DOX to DOXol is predominantly catalyzed in the cytosol [5]. Kinetic data for reduction of DOX by purified CBR1 and AKR1A1 showed that CBR1 had lower Km and higher turnover number (kcat), and consequently higher catalytic efficiency (kcat/Km), with respect to AKR1A1 ([14], Table 1). This can be explained in terms of the shorter distance
152
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
between the substrate carbonyl group in C13 and the enzyme active site as shown by molecular docking and dynamics simulations and reported in Figs. 3e7. As the hydroxyl group in C14 is replaced by a hydrogen (DAUN) the anthracycline is more tightly bound to the enzymeecoenzyme complex (Figs. 6 and 7) and this is expected to increase the catalytic activity of AKR1A1 which is consistent with the experimental observation that the specific activity of AKR1A1 toward DAUN was higher than that toward DOX ([18], Table 1). It is also interesting to compare the DAUN/CBR1 and DOX/CBR1 binding modeled in the present study with two different experimental measurements which support that DAUN is a better substrate than DOX for CBR1 on the basis of significantly higher kcat/Km values ([13,15], Table 1). Our computational results point out the role of a hydrophobic group in C14 in imparting stability to the productive enzymeesubstrate complex. This stabilization would be impaired for DOX which in the course of the simulation interacts less efficiently with the residues of the catalytic site. In summary, our modeled binding modes of daunorubicin and doxorubicin with cytosolic reductases can reasonably explain all of the previously reported experimental data. 3. Conclusions Anticancer DAUN and DOX reductions are associated with increasing risk of chronic cardiotoxicity. One of the factors underlying the interpatient variation in cardiotoxicity may be linked to cytosolic reductases responsible for metabolizing DOX and DAUN. Several reductase enzymes are involved in the metabolism of these drugs, and distinguishing their individual contributions to the total metabolism is a challenging task. It is with these considerations that we present our molecular modeling findings focusing on the interaction of DOX and DAUN with CBR1 and AKR1A1, the most efficient cytosolic reductases acting on anthracycline antibiotics. We describe a combined study in which homology modeling, docking and MD simulations were employed in order to rationalize the question of how the two different substrates influence the selectivity of CBR1 and AKR1A1. Kinetic studies in vitro indicate that CBR1 and AKR1A1 are less active in metabolizing DOX, with respect to DAUN, by the decrease of kcat/Km, the best kinetic parameter to assess overall substrate efficiency, and, in all cases, CBR1 results a more efficient enzyme (Table 1). The ligandeenzyme interaction models, chosen on the basis of both catalytically productive orientation and docking energy score and further evaluated for dynamic stability of the complexes throughout the MD simulations, are fully consistent with the commonly accepted mechanism of catalysis of CBR1 and AKR1A1 and with previously reported kinetic studies. We also show that the more hydrophobic daunorubicin side chain in C9 may have an important role in the substrate recognition mechanism of these enzymes and that a hydrophobic group in C14 confers stability to the enzymeesubstrate complex. Given that inhibitors of the intramyocardial formation of DOXol and DAUNol metabolites might prove useful to mitigate cardiotoxicity, our results have therefore implications for drug discovery. In future outlooks, the binding models presented here, might serve to design a minimal pharmacophore for the specific enzyme-substrate recognition with the aim to design a novel class of inhibitors of anthracycline reductases.
PP (3-(1-tert-butyl-4-amino-1H-pyrazolo[3,4-d]pyrimidin-3-yl) phenol) molecule. To generate docking targets, inhibitors and crystallographic water molecules were removed from PDB structure. The nicotinamide mononucleotide unit of the cofactor was reduced at C4. The crystallographic structure of AKR1A1 (PDB ID: 2ALR [30]) is missing residues 214e225 and a 3D model structure of AKR1A1 was therefore built by homology modeling. AKR1A1 displays high sequence identity (94%) to porcine aldehyde reductase. Thus, the crystal structure of porcine aldehyde reductase in complex with fidarestat (PDB ID 3H4G [31]) was used as structural template. The sequence alignment produced by ClustalW was used, in conjunction with the program MODELLER [32] implemented in DiscoveryStudio v3.0 (Accelrys Inc.), to generate a set of 50 models of AKR1A1 in complex with fidarestat. The homology model was compared to the 2ALR crystal structure by calculating the rootmean square deviation (RMSD) of the a-carbons in DiscoveryStudio. Only the residues present in the crystal structure were compared to the homology model. RMSD values were also calculated for heavy atoms within 3.5 Å of the cocrystallized fidarestat (residues 21, 48, 49, 112, 113, 298, 299, 300). The stereochemical accuracy of the model was also checked by means of PROCHECK v.3.5.4 [33] and VERIFY3D [34] software. DOX and DAUN structures were extracted from their DNA complexes (PDB ID: 151D [35] and 1D33 [36], respectively) and subjected to 1000 steps of steepest descent energy minimization in vacuo, using GROMACS v. 4.5 [37] and the GROMOS96 43a1 force field [38]. The protonation state 30 -NH3 þ in the daunosamine ring was used. EPI and MEN 10755 were built starting from the structure of DOX using DiscoveryStudio v3.0 (Accelrys Inc.). 4.2. Docking calculations Flexible-ligand docking was carried out using AutoDock Vina 1.1 [39]. This software has been shown to produce, together with an increased efficiency in predicting the experimental binding poses and energies, a 2 orders of magnitude speed-up compared with AutoDock 4 [40]. Vina uses as input files the 3D coordinates of both ligand and receptor, which must be converted into the appropriate format using the AutoDockTools (ADT) program [41]. Polar hydrogens were added by using the Hydrogens module in ADT; after that, Gasteiger partial charges were assigned [42]. For all docking calculations, a grid box size of 30 30 30 points (spacing between the grid points of 1 Å) was used, centered on the mass center of the crystallographic ligand and encompassing all active site atoms. The AutoDock Vina parameter “Exhaustiveness”, which determines how comprehensively the program searches for the lowest energy was set to the value of 500. All other parameters were used as defaults. All rotatable bonds within the ligand were allowed to rotate freely, and the receptor was considered rigid. All computations were carried out on a HP XW8600 workstation running Red Hat Enterprise Linux 5. The resulting docking poses were analyzed with DiscoveryStudio to detect catalytically productive orientation of the substrate within the binding site. Schematic two-dimensional representations of the docking results were produced using LIGPLOT [43]. 4.3. Validation of the docking protocol
4. Experimental section 4.1. Coordinates of proteins and substrates Structure of Protein Data Bank (PDB) entry 1WMA [26] was used for docking simulations of CBR1. The 1WMA structure contains one molecule of CBR1 complexed with an NADPþ and with a hydroxy-
The structures of CBR1 bound to hydroxy-PP (PDB ID: 1WMA) and porcine aldehyde reductase (PDB ID: 3H4G) bound to fidarestat were used for docking validation (correctly identifying the binding site and scoring the receptoreligand interactions). Hydroxy-PP and fidarestat were docked to the targets containing CBR1eNADPH and porcine aldehyde reductase-NADPH, respectively. Ligand docked
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
poses were in close agreement with the crystallographically determined positions. The RMSD of the heavy atoms was only 1.45 Å for hydroxy-PP in CBR1 and 1.12 Å for fidarestat in porcine aldehyde reductase. These control simulations demonstrate that AutoDock Vina is able to generate docking models in the chosen crystal structures and can be used to predict DOX and DAUN binding modes in cytosolic reductases. 4.4. Molecular dynamics calculations Enzymeesubstrate complexes and free enzyme dynamics were conducted using the same method. Simulations and analysis of the trajectories were carried out using GROMACS 4.5 [37] and the GROMOS96 43a1 force field [38]. The topology of the protein and NADPH was created using GROMACS utilities, whereas the topology of the ligand was created using the PRODRG server [44], available at http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg. All starting structures were immersed in a triclinic water box under periodic boundary conditions using a 1.0 nm minimum distance from the protein to the box faces. Each system was then neutralized by Cl and Naþ counterions (CBR1 and AKR1A1 systems, respectively) that were added at random positions to the bulk solvent. The final systems consisted of 2690 and 3258 protein atoms (CBR1 and AKR1A1, respectively) surrounded by w11,000 water molecules. Following steepest descents minimization, the systems were equilibrated under NVT conditions for 100 ps at 300 K, followed by 100 ps under NPT conditions, applying position restraints to the protein and ligand atoms. Finally, all restraints were removed and the molecular dynamics were run for 20 ns at 300 K. The temperature was maintained close to its reference value (300 K) by applying the V-rescale [45] thermostat with a coupling constant of 0.1 ps whereas for the constant pressure (1 atm, isotropic coordinate scaling) the ParrinelloeRahman barostat [46] was used with a relaxation time of 2.0 ps. van der Waals interactions were modeled using 6e12 Lennard-Jones potentials with a 1.4 nm cutoff. Long-range electrostatic interactions were calculated using Particle Mesh Ewald method, with a cutoff for the real space term of 0.9 nm. Covalent bonds were constrained using the LINCS algorithm. The time step employed was 2 fs and the coordinates were saved every 2 ps for analysis which was carried out using the standard GROMACS tools. Acknowledgments Computational resources for MD simulations were provided by CASPUR e Inter-University Consortium for Supercomputing application for University and Research (Standard HPC Grant). References [1] R.B. Weiss, The anthracyclines: will we ever find a better doxorubicin? Semin. Oncol. 19 (1992) 670e686. [2] R.C. Young, R.F. Ozols, C.E. Myers, The anthracycline antineoplastic drugs, N. Engl. J. Med. 305 (1981) 139e153. [3] S.B. Kaye, The multidrug resistance phenotype, Br. J. Cancer 58 (1988) 691e694. [4] L.J. Steinherz, Anthracycline-induced cardiotoxicity, Ann. Intern. Med. 126 (1997) 827. author reply 827e828. [5] G. Minotti, P. Menna, E. Salvatorelli, G. Cairo, L. Gianni, Anthracyclines: molecular advances and pharmacologic developments in antitumor activity and cardiotoxicity, Pharmacol. Rev. 56 (2004) 185e229. [6] R. Zucchi, R. Danesi, Cardiac toxicity of antineoplastic anthracyclines, Curr. Med. Chem. Anticancer Agents 3 (2003) 151e171. [7] R.D. Olson, P.S. Mushlin, D.E. Brenner, S. Fleischer, B.J. Cusack, B.K. Chang, et al., Doxorubicin cardiotoxicity may be caused by its metabolite, doxorubicinol, Proc. Natl. Acad. Sci. U. S. A. 85 (1988) 3585e3589. [8] B.J. Cusack, P.S. Mushlin, L.D. Voulelis, X. Li, R.J. Boucek Jr., R.D. Olson, Daunorubicin-induced cardiac injury in the rabbit: a role for daunorubicinol? Toxicol. Appl. Pharmacol. 118 (1993) 177e185. [9] R.L. Felsted, N.R. Bachur, Mammalian carbonyl reductases, Drug Metab. Rev. 11 (1980) 1e60.
153
[10] A. Mordente, G. Minotti, G.E. Martorana, A. Silvestrini, B. Giardina, E. Meucci, Anthracycline secondary alcohol metabolite formation in human or rabbit heart: biochemical aspects and pharmacologic implications, Biochem. Pharmacol. 66 (2003) 989e998. [11] T. Matsunaga, S. Shintani, A. Hara, Multiplicity of mammalian reductases for xenobiotic carbonyl compounds, Drug Metab. Pharmacokinet. 21 (2006) 1e18. [12] H. Ohara, Y. Miyabe, Y. Deyashiki, K. Matsuura, A. Hara, Reduction of drug ketones by dihydrodiol dehydrogenases, carbonyl reductase and aldehyde reductase of human liver, Biochem. Pharmacol. 50 (1995) 221e227. [13] A. Slupe, B. Williams, C. Larson, L.M. Lee, T. Primbs, A.J. Bruesch, et al., Reduction of 13-deoxydoxorubicin and daunorubicinol anthraquinones by human carbonyl reductase, Cardiovasc. Toxicol. 5 (2005) 365e376. [14] N. Kassner, K. Huse, H.-J. Martin, U. Gödtel-Armbrust, A. Metzger, I. Meineke, et al., Carbonyl reductase 1 is a predominant doxorubicin reductase in the human liver, Drug Metab. Dispos. 36 (2008) 2113e2120. [15] O.S. Bains, M.J. Karkling, T.A. Grigliatti, R.E. Reid, K.W. Riggs, Two nonsynonymous single nucleotide polymorphisms of human carbonyl reductase 1 demonstrate reduced in vitro metabolism of daunorubicin and doxorubicin, Drug Metab. Dispos. 37 (2009) 1107e1114. [16] J.G. Blanco, W.M. Leisenring, V.M. Gonzalez-Covarrubias, T.I. Kawashima, S.M. Davies, M.V. Relling, et al., Genetic polymorphisms in the carbonyl reductase 3 gene CBR3 and the NAD(P)H:quinone oxidoreductase 1 gene NQO1 in patients who developed anthracycline-related congestive heart failure after childhood cancer, Cancer 112 (2008) 2789e2795. [17] R.H. Takahashi, O.S. Bains, T.A. Pfeifer, T.A. Grigliatti, R.E. Reid, K.W. Riggs, Aldo-keto reductase 1C2 fails to metabolize doxorubicin and daunorubicin in vitro, Drug Metab. Dispos. 236 (2008) 991e994. [18] T. O’connor, L.S. Ireland, D.J. Harrison, J.D. Hayes, Major differences exist in the function and tissue-specific expression of human aflatoxin B1 aldehyde reductase and the principal human aldo-keto reductase AKR1 family members, Biochem. J. 343 (1999) 487e504. [19] H.-J. Martin, U. Breyer-Pfaff, V. Wsol, S. Venz, S. Block, E. Maser, Purification and characterization of akr1b10 from human liver: role in carbonyl reduction of xenobiotics, Drug Metab. Dispos. 34 (2006) 464e470. [20] A. Mordente, E. Meucci, A. Silvestrini, G.E. Martorana, B. Giardina, New developments in anthracycline-induced cardiotoxicity, Curr. Med. Chem. 16 (2009) 1656e1672. [21] G. Minotti, S. Licata, A. Saponiero, P. Menna, A.M. Calafiore, G. Di Giammarco, et al., Anthracycline metabolism and toxicity in human myocardium: comparisons between doxorubicin, epirubicin, and a novel disaccharide analogue with a reduced level of formation and [4Fe-4S] reactivity of its secondary alcohol metabolite, Chem. Res. Toxicol. 13 (2000) 1336e1341. [22] W. Fiedler, N. Tchen, J. Bloch, P. Fargeot, R. Sorio, J.B. Vermorken, et al., A study from the EORTC new drug development group: open label phase II study of sabarubicin (MEN-10755) in patients with progressive hormone refractory prostate cancer, Eur. J. Cancer 42 (2006) 200e204. [23] E. Salvatorelli, P. Menna, M. Lusini, E. Covino, G. Minotti, Doxorubicinolone formation and efflux: a salvage pathway against epirubicin accumulation in human heart, J. Pharmacol. Exp. Ther. 329 (2009) 175e184. [24] G. Minotti, S. Recalcati, A. Mordente, G. Liberi, A.M. Calafiore, C. Mancuso, P. Preziosi, G. Cairo, The secondary alcohol metabolite of doxorubicin irreversibly inactivates aconitase/iron regulatory protein-1 in cytosolic fractions from human myocardium, FASEB J. 12 (1998) 541e552. [25] V. Gonzalez-Covarrubias, J.L. Kalabus, J.G. Blanco, Inhibition of polymorphic human carbonyl reductase 1 (CBR1) by the cardioprotectant flavonoid 7monohydroxyethyl rutoside (monoHER), Pharm. Res. 25 (2008) 1730e1734. [26] M. Tanaka, R. Bateman, D. Rauh, E. Vaisberg, S. Ramachandani, C. Zhang, et al., An unbiased cell morphology-based screen for new, biologically active small molecules, PLoS Biol. 3 (2005) e128. [27] U. Oppermann, C. Filling, M. Hult, N. Shafqat, X. Wu, M. Lindh, et al., Shortchain dehydrogenases/reductases (SDR): the 2002 update, Chem. Biol. Interact 143e144 (2003) 247e253. [28] J.M. Jez, M.J. Bennett, B.P. Schlegel, M. Lewis, T.M. Penning, Comparative anatomy of the aldo-keto reductase superfamily, Biochem. J. 326 (Pt 3) (1997) 625e636. [29] B.P. Schlegel, J.M. Jez, T.M. Penning, Mutagenesis of 3 alpha-hydroxysteroid dehydrogenase reveals a «push-pull» mechanism for proton transfer in aldo-keto reductases, Biochemistry 37 (1998) 3538e3548. [30] O. El-Kabbani, N.C. Green, G. Lin, M. Carson, S.V. Narayana, K.M. Moore, et al., Structures of human and porcine aldehyde reductase: an enzyme implicated in diabetic complications, Acta Crystallogr. D Biol. Crystallogr. 50 (1994) 859e868. [31] O. El-Kabbani, V. Carbone, C. Darmanin, M. Oka, A. Mitschler, A. Podjarny, et al., Structure of aldehyde reductase holoenzyme in complex with the potent aldose reductase inhibitor fidarestat: implications for inhibitor binding and selectivity, J. Med. Chem. 48 (2005) 5536e5542. [32] A. Sali, T.L. Blundell, Comparative protein modelling by satisfaction of spatial restraints, J. Mol. Biol. 234 (1993) 779e815. [33] A.L. Morris, M.W. MacArthur, E.G. Hutchinson, J.M. Thornton, Stereochemical quality of protein structure coordinates, Proteins 12 (1992) 345e364. [34] R. Lüthy, J.U. Bowie, D. Eisenberg, Assessment of protein models with threedimensional profiles, Nature 356 (1992) 83e85. [35] L.A. Lipscomb, M.E. Peek, F.X. Zhou, J.A. Bertrand, D. VanDerveer, L.D. Williams, Water ring structure at DNA interfaces: hydration and dynamics of DNAanthracycline complexes, Biochemistry 33 (1994) 3649e3659.
154
D. Pirolli et al. / European Journal of Medicinal Chemistry 56 (2012) 145e154
[36] A.H. Wang, Y.G. Gao, Y.C. Liaw, Y.K. Li, Formaldehyde cross-links daunorubicin and DNA efficiently: HPLC and X-ray diffraction studies, Biochemistry 30 (1991) 3812e3815. [37] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, GROMACS 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theory Comput. 4 (2008) 435e447. [38] W.R.P. Scott, P.H. Hünenberger, I.G. Tironi, A.E. Mark, S.R. Billeter, J. Fennen, et al., The GROMOS biomolecular simulation program package, J. Phys. Chem. A 103 (1999) 3596e3607. [39] O. Trott, A.J. Olson, AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading, J. Comput. Chem. 31 (2010) 455e461. [40] R. Huey, G.M. Morris, A.J. Olson, D.S. Goodsell, A semiempirical free energy force field with charge-based desolvation, J. Comput. Chem. 28 (2007) 1145e1152.
[41] M.F. Sanner, Python: a programming language for software integration and development, J. Mol. Graph. Model. 17 (1999) 57e61. [42] J. Gasteiger, Iterative partial equalization of orbital electronegativity-a rapid access to atomic charges, Tetrahedron 36 (1980) 3219e3228. [43] A.C. Wallace, R.A. Laskowski, J.M. Thornton, LIGPLOT e a program to generate schematic diagrams of protein ligand interactions, Protein Eng. 8 (1995) 127e134. [44] A.W. Schüttelkopf, D.M.F. van Aalten, PRODRG: a tool for high-throughput crystallography of protein-ligand complexes, Acta Crystallogr. D Biol. Crystallogr. 60 (2204) 1355e1363. [45] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (2007) 014101. [46] M. Parrinello, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182.