Chem Biol Drug Des 2015 Research Article
Characterizing the Free-Energy Landscape of MDM2 Protein–Ligand Interactions by Steered Molecular Dynamics Simulations Guodong Hu1,2,*, Shicai Xu1,2 and Jihua Wang1,2,* 1
Shandong Provincial Key Laboratory of Functional Macromolecular Biophysics, Dezhou University, Dezhou, China 2 College of Physics and Electronic Information, Dezhou University, Dezhou, China *Corresponding authors: Guodong Hu, [email protected]; Jihua Wang, [email protected] Inhibition of p53–MDM2 interaction by small molecules is considered to be a promising approach to re-activate wild-type p53 for tumor suppression. Several inhibitors of the MDM2–p53 interaction were designed and studied by the experimental methods and the molecular dynamics simulation. However, the unbinding mechanism was still unclear. The steered molecular dynamics simulations combined with Brownian dynamics fluctuation–dissipation theorem were employed to obtain the free-energy landscape of unbinding between MDM2 and their four ligands. It was shown that compounds 4 and 8 dissociate faster than compounds 5 and 7. The absolute binding free energies for these four ligands are in close agreement with experimental results. The open movement of helix II and helix IV in the MDM2 protein-binding pocket upon unbinding is also consistent with experimental MDM2unbound conformation. We further found that different binding mechanisms among different ligands are associated with H-bond with Lys51 and Glu25. These mechanistic results may be useful for improving ligand design. Key words: MDM2 protein, steered molecular dynamics simulations, unbinding mechanism Received 12 April 2015, revised 20 May 2015 and accepted for publication 25 May 2015
The human version of mouse double minute protein 2 (MDM2) inhibits tumor suppressor activity of p53 through binding of its N-terminal domain to the transactivation domain of p53 (1–4). Disruption of the binding between MDM2 and p53 would enhance the ability of p53 to suppress tumor growth. As a result, MDM2 is an important anticancer target for small-molecule inhibitor design (5–8). ª 2015 John Wiley & Sons A/S. doi: 10.1111/cbdd.12598
Using a structure-based approach, Ding et al. (9) successfully designed several inhibitors of the MDM2–p53 interaction. The most potent ligand 8 has a Ki value of 3 nM binding to MDM2 and is 2000 times more potent than the natural p53 peptide (residues 13–19). These inhibitors were designed to mimic the amino acid residues in the short a-helical region of p53 (Phe19, Trp23, and Leu26) that interacts with MDM2. Molecular dynamics simulations coupled with free-energy calculations by molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) confirmed those ligands as residue mimicries although not all ligands mimicked the particular residue as designed, such as ligands 4 and 8 with a 2-morpholin-4yl-ethylamine group. However, the information about the unbinding of ligands from the binding pocket of MDM2 was still unclear (10). The purpose of this study was to examine the energy landscape of ligand–MDM2 interactions by steered molecular dynamics (SMD) simulations. SMD takes inspiration from single-molecule pulling experiments (11) and dissociates a complex structure by a pulling force (12,13). The non-equilibrium dynamics of the system under a pulling force can extract the free-energy difference between the bound state and the unbound state by measuring the work along the transition paths and map out the freeenergy landscape in terms of the potential of mean force (PMF) (14,15) with high precision and efficiency (16–18). SMD simulations have become widely used in studying many biochemical processes including the transportation of ions and organic molecules across membrane channels (16,19–22) and unfolding/folding mechanism of proteins (23,24) as well as the mechanisms of protein–ligand binding (12,17,25). In this study, we combined SMD with Brownian dynamics fluctuation–dissipation theorem (BD-FDT) to calculate the interaction energy landscape of MDM2 with four known MDM2 ligands (compounds 4, 5, 7, and 8, Figure 1) (9). The largest difference between ligand 8 and ligand 5 is the 2-morpholin-4yl-ethylamine group. The van der Waals energies, electrostatic energies, and hydrogen bonds were analyzed in detail along the pulling paths of SMD simulation. The binding affinities of these four ligands from energy landscape studies are in close agreement with experimental measurements. 1
Hu et al.
Figure 1: Molecular structures of four ligands which are labeled with their name as in reference (7) (compounds 4, 5, 7 and 8).
Model and Method Systems setup The starting structures for MDM2–ligand complexes without water molecules and counter ions were obtained from the last conformation of 3-nseconds equilibrated MD simulations of our previous work (10). To minimize steric
crashes during SMD simulations, all initial complex structures were rotated so that the orifice of MDM2’s binding pocket points to +z direction (Figure 2). The complexes were then solvated in a rectangular water box using VMD software (26). The minimum distances from the surface of complexes to the edges of boxes are 15 A in x, y, and –z directions and 25 A in the +z direction. Afterward, an
Figure 2: The orientations of ligands in the binding pocket of MDM2 from +z direction to the complex. The ligands are shown in stick and ball representations; the pulled atoms are shown in large balls.
2
Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes
appropriate number of counterions were added to neutralize the charge of systems. They were placed at the grids with the largest positive coulombic potentials around the complexes using the Tleap module implemented in a AMBER10 package . The standard AMBER force field for bioorganic systems (FF03) (27) was utilized to describe the proteins and water molecule (TIP3P). The parameters for ligands were obtained previously followed a standard AMBER procedure for force-field parameter and charge assignment through quantum calculations (10).
Molecular dynamics (MD) simulation and steered molecular dynamics (SMD) simulation The MD simulations and SMD simulations were performed using NAMD 2.8 package (15). We employed periodic boundary conditions in all three dimensions and utilized particle mesh Ewald (PME) method for treating long-range electrostatic interactions. The cutoff distance for longrange interactions was set to 12 A with a switching distance of 10 A. The time step was 2 fseconds for shortrange and 4 fseconds for long-range interactions. The covalent bonds of all hydrogen atoms were fixed to their equilibrium length by SHAKE algorithm (28). The temperature is maintained at 300 K using Langevin thermostat (29) with the damping constant chosen as 5/pseconds and pressure at 1 bar. To remove possible steric crashes between water molecules and complex structures, the energy of each system was minimized for 2000 steps, during which locations of protein and ligand atoms were restrained by a harmonic potential with a force constant of 0.2 kcal/(mol 9 A2). Afterward, 5-nseconds isothermal isobaric ensemble (NPT)-MD simulations were performed without any restraints. In our SMD, a ligand was pulled from a bound state to a dissociated state. The mass center of a ligand was initially set at the origin of co-ordinates. To maintain the location of MDM2 during the SMD simulations, we restrained the z-coordinates of several backbone carbon atoms of MDM2 that are far away from the binding pocket (Figure 3A). Four MDM2–ligand complex systems were equilibrated for 150 000 steps for each segment without pulling forces and for 50 000 steps with pulling forces.
Potential of mean force (PMF) from steered molecular dynamics (SMD) simulation Because the four ligands studied in this study have more than seventy atoms each, one center for a ligand is unlikely to sample all relevant events. As a result, each ligand is represented by two centers (Figure 2). Both centers were steered along z direction. The pulling speed was set at 0.01 A/pseconds at the z direction, and a ligand can move freely in x and y directions. To reduce the impact of pulling on MDM2, we divided the pulling path into 16 segChem Biol Drug Des 2015
A
B
Figure 3: (A) Pulling compound 8 from its bound state to an unbound state. MDM2 (a new cartoon and a surface representation), compound 8 (a ball representation), and the pulling path (the red line) are shown along with the displacements along z-axis and in xy plane as labeled. The restrained atoms are shown in a ball representation with gray color. (B) Potential of mean force for four ligands (compounds 4, 5, 7, and 8) as a function of the displacement from their binding site along the pulling path when two centers were steered away from the MDM2 protein.
ments along the z direction with 1 A for each segment. For each segment, we performed two types of SMD simulations: one for pulling away (denoted as forward) and one for pulling back to (denoted as reverse) the binding site. For each segment, two end states are denoted as A and B, respectively. Eight forward and reverse pulling paths between A and B states are sampled during which the work done to the system was recorded. The Gibbs freeenergy difference between state A and an intermediate state r (namely, the PMF) was given by the Brownian dynamics fluctuation–dissipation theorem (BD-FDT) (14,30) as follows: hexp½WA!r =2kB TiF DG ¼ kB T ln hexp½Wr!A =2kB TiR
where kB and T are the Boltzmann constant and the absolute temperature, respectively. The brackets with subscript F/R represent the statistical average over the forward/ reverse paths, respectively. WA?r is the work done to the 3
Hu et al.
Figure 4: Work done along the forward and reverse pulling paths. In the forward direction, the atoms of a ligand were pulled away from the binding site. In the reverse direction, the atoms of a ligand were pulled back to the binding site.
system along a forward path when the ligand is steered from A to r. Wr?A (=WB?A – WB?r) is the work done to the system along a reverse path when the ligand is steered from r to A. Data analysis We define that a ligand and MDM2 are in contact if an atom of MDM2 is within 7.0 A of any atoms of a ligand. The hydrogen bonds (H-bonds) are defined by the distance between acceptor and donor atoms ( 13 A, indicating that a 25- A water box for the +z direction is sufficient for the SMD simulations. Compounds 4 and 8, however, reached plateau at about 10 A, compared to 12 A for compounds 5 and 7; that is, the former can dissociate more quickly than the latter. Compound 5 reveals an interesting intermediate state at around 4–6 A. On the other hand, the PMF for compound 7 is made of two phases: a rapid increase phase from 0 to 3 A and a steady but slower increase phase after 3 A. Thus, the binding thermodynamics of four ligands are different. To understand different mechanisms for different compounds observed above, we calculated the interacting Table 1: The binding free energies for four complexesa Compound b
Figure 4 compares the work required for forward and reverse directions. The figure shows that the pulling process is largely irreversible on the basis of the following details. WA?z is not equal to Wz?A in most of the segments, especially in the first few segments of pulling. WA? z would be equal to Wz?A if the pulling was reversible. 4
MM-PBSA BD-FDT (this work) Expc
4
5
7
8
16.00 9.93 10.83
10.74 9.53 10.19
12.36 9.59 10.1
19.64 10.88 11.7
a
All values are give in kcal/mol. The binding free energies were calculated using the MM-PBSA method (8). c The experimental binding free energies were calculated using Ki (7). b
Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes A
B
C
Figure 5: The averaged energies as a function of ligand displacement along z-axis were (A) for van der Waals energies and (B) for electrostatic energies. (C) The averaged number of Hbonds between ligands and MDM2 as a function of ligand displacement along z-axis; the key values were labeled with ligand name and letters.
energy terms for the structures along the pulling path. For a given z value, van der Waals energies and the electrostatic energies were averaged over 16 structures extracted from eight forward SMD trajectories and eight reverse SMD trajectories. Figure 5A shows that the van der Waals energies increase smoothly as a function of the ligand displacement along the z-axis for four complexes. The magnitudes of van der Waals energies for four complexes at z = 0 have the same rank order as their experimental binding affinity values. Smooth increase of van der Waals energies as z increases indicates lack of steric hindrance for all ligands transiting from a abound state to an unbound state. However, large fluctuations of electrostatic energies are observed (Figure 5B). The minima for electrostatic Chem Biol Drug Des 2015
energies are located between 5 A to 8 A for compound 4, 8 A to 12 A for compound 5, and 8 A to 10 A for compound 7. To have a better understanding of large fluctuation in electrostatic energies, we calculated the number of H-bonds between a ligand and MDM2 (Figure 5C). Although the average number of protein–ligand H-bonds is mostly smaller than one, all the large dips in electrostatic energies for compounds 4, 5, and 7 correspond to a peak in the number of H-bonds. This indicates the importance of hydrogen bonding energies in overall electrostatic energies, despite that the number of H-bonds is mostly 0.5) with MDM2 until z > 8 A. The H-bonds were formed between its N3 H6 and the oxygen atoms of side chain of the flexible N-terminal residue Glu25 which allows a longer distance from compound 5 to MDM2 (Figure 7C). These H-bonds are semistable as their probability is about 0.5 (Figure 5C).
Unbinding mechanisms of compounds 7 and 8 Unlike compounds 4 and 5, the H-bonds between compounds 7/8 and MDM2 decay more slowly. The number of H-bond between MDM2 and compound 7 decreases smoothly from z = 0 to z = 6 A and then increases with a small peak at z = 8 A. For compound 8, its H-bond to MDM2 starts to decrease only at z = 6 A and becomes nearly zero at z > 8 A. Compounds 7 and 8 both had an H-bond to the oxygen atom of backbone of Leu54 (Figures 8B and 9) as with compounds 4 and 5. However, unlike compounds 4 and 5, geometry of compounds 7 and 8 allowed a large displacement as shown in Figure 9. Unlike compound 8, the atoms N3 H6 of ligand 7 can form a new H-bond with the oxygen atom of backbone of Lys51 at z = 8 A. Figure 9 shows the overlapped four structures extracted from the SMD simulation with the center of ligand at the position shown in Figure 5C. In summary, the H-bonds between nitrogen atom of a ligand and the backbone oxygen atom of Leu54 played 5
Hu et al.
B A
D
C
Figure 6: Schematic view of the key distances between ligand 4 and MDM2. MDM2 was shown in new cartoon representation and colored by secondary structure. Ligand was shown in stick representation and colored by name of atoms. The key atoms were shown in ball representation. The distance was labeled with dash line. The plots A, B, C, and D were structures for the center of ligand 4 at 0 A, 3.93 A, 5.36 A, and 7.11 A, respectively.
A
B
D
a key role for MDM2–ligand binding, consistent with previous work (9,10). Their stabilities are different for different compounds. After breaking the H-bond to Leu54, compounds 4 and 7 can form H-bond with Lys51, 6
C
Figure 7: Schematic view was similar to Figure 6 except for ligand 5. The plots A, B, C, and D were structures for the center of ligand 5 at 0 A, 3.06 A, 10.63 A, and 12.3 A, respectively.
whereas compound 5 forms H-bond with Glu25, but not compound 8. Thus, Lys51 and Glu25 are important for ligand dissociations from MDM2 for three of four ligands studied. Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes
Figure 8: Schematic view was similar to Figure 6 except for ligand 7. The plots A, B, C, and D were structures for the center of ligand 7 at 0 A, 3.74 A, 6.56 A, and 8.25 A, respectively.
A
B
C
D
Figure 9: Four structures of 8a (0 A), 8b (4.69 A), 8c (7.55 A), and 8d (8.54 A) were overlapped according to MDM2. The key section was enlarged and labeled.
The conformation of the binding cleft To study the effect of ligand dissociations on the dynamics and conformations in the binding pocket of the MDM2, RMSFs versus the residue number for MDM2 in the bound (z ≤ 1 A) and unbound states (15 A < z ≤ 16 A) were analyzed (Figure 10). The RMSFs of the most residues in MDM2 in the unbound state were larger than those in the bound state, the residues 41–56, and residues 91–100 that consist of loop regions and portions of helix II and helix IV. We calculated the averaged distances of four compounds between the center of helix II and helix IV in the bound and unbound states. The averaged fluctuation in the unbound state increases by about 0.43 A than the bound state, consistent with 0.5 A and 0.7 A increase upon p53 unbinding (31–33). Chem Biol Drug Des 2015
Conclusions Steered molecular dynamics combined with free-energy calculations is emerging as a promising tool in structuredbased drug design. In this study, SMD combined with BDFDT was used to study the interaction of four ligands with the MDM2 protein. The PMFs of these four ligands indicate that compounds 4 and 8 can reach an unbound state faster than compounds 5 and 7. This is likely because compounds 5 and 7 can have H-bonds with Lys51 at a large displacement from MDM2. The binding affinities from the PMF were in good agreement with the experimental data. This work provides the atomistic details of the unbinding mechanism, revealing transient interactions at accessory binding sites which can be exploited in-depth for inhibitor design for MDM2. 7
Hu et al.
Figure 10: Root mean square fluctuation of backbone atoms versus residue number of MDM2. The regions of four helices of MDM2 were labeled.
Acknowledgments This work was partially supported by funding from the National Natural Science Foundation of China (61271378 and 11447004), the Natural Science Foundation of Shandong Province (ZR2014JL006 and ZR2012CL09), and A Project of Shandong Province Higher Educational Science and Technology Program (J14LJ05). We would like to thank Yaoqi Zhou for his scientific support.
Conflict of Interests The authors declare that they have no conflict of interest.
References 1. Vogelstein B., Lane D., Levine A.J. (2000) Surfing the p53 network. Nature;408:307–310. 2. Levine A.J. (1997) p53, the cellular gatekeeper for growth and division. Cell;88:323–331. 3. Zhang X.P., Liu F., Cheng Z., Wang W. (2009) Cell fate decision mediated by p53 pulses. Proc Natl Acad Sci U S A;106:12245–12250. 4. Zhang X.P., Liu F., Wang W. (2011) Two-phase dynamics of p53 in the DNA damage response. Proc Natl Acad Sci U S A;108:8990–8995. 5. Vassilev L.T., Vu B.T., Graves B., Carvajal D., Podlaski F., Filipovic Z., Kong N., Kammlott U., Lukacs C., Klein C. (2004) In vivo activation of the p53 pathway by small-molecule antagonists of MDM2. Science;303:844–848. 6. Ding K., Lu Y., Nikolovska-Coleska Z., Qiu S., Ding Y., Gao W., Stuckey J., Krajewski K., Roller P.P., Tomita Y., Parrish D.A., Deschamps J.R., Wang S. (2005) Structure-based design of potent non-peptide MDM2 inhibitors. J Am Chem Soc;127:10130–10131. 7. Vassilev L.T. (2005) p53 activation by small molecules: application in oncology. J Med Chem;48:4491–4499. 8
8. Hang Y., Gui-in L., Hyung Soon P., Gregory A.P., Johanna M.R., Said M.S., Andrew D.H. (2005) Terphenyl-based helical mimetics that disrupt the p53/HDM2 interaction13. Angew Chem Int Ed;44:2704–2707. 9. Ding K., Lu Y., Nikolovska-Coleska Z., Wang G., Qiu S., Shangary S., Gao W., Qin D., Stuckey J., Krajewski K., Roller P.P., Wang S. (2006) Structure-based design of spiro-oxindoles as potent, specific small-molecule inhibitors of the MDM2-p53 interaction. J Med Chem;49:3432–3435. 10. Hu G., Wang D., Liu X., Zhang Q. (2010) A computational analysis of the binding model of MDM2 with inhibitors. J Comput Aided Mol Des;24:687–697. 11. Grubmuller H., Heymann B., Tavan P. (1996) Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force. Science;271:997–999. 12. Patel J.S., Berteotti A., Ronsisvalle S., Rocchia W., Cavalli A. (2014) Steered molecular dynamics simulations for studying protein-ligand interaction in cyclindependent kinase 5. J Chem Inf Model;54:470–480. 13. Li M.S., Mai B.K. (2012) Steered molecular dynamics – a promising tool for drug design. Curr Bioinform;7:342. 14. Chen L.Y. (2008) Nonequilibrium fluctuation-dissipation theorem of Brownian dynamics. J Chem Phys;129:144113. 15. Allen T.W., Andersen O.S., Roux B. (2006) Molecular dynamics – potential of mean force calculations as a tool for understanding ion permeation and selectivity in narrow channels. Biophys Chem;124:251–267. 16. Chen L.Y. (2013) Glycerol modulates water permeation through Escherichia coli aquaglyceroporin GlpF. Biochim Biophys Acta;1828:1786–1793. 17. Chen L.Y. (2011) Exploring the free-energy landscapes of biological systems with steered molecular dynamics. Phys Chem Chem Phys;13:6176–6183. 18. Hu G., Chen L.Y. (2010) In silico experiments of single-chain antibody fragment against drugs of abuse. Biophys Chem;153:97–103. 19. Giorgino T., De Fabritiis G. (2011) A high-throughput steered molecular dynamics study on the free energy profile of ion permeation through gramicidin A. J Chem Theory Comput;7:1943. 20. Henin J., Tajkhorshid E., Schulten K., Chipot C. (2008) Diffusion of glycerol through Escherichia coli aquaglyceroporin GlpF. Biophys J;94:832–839. 21. Mai B.K., Viet M.H., Li M.S. (2010) Top leads for swine influenza A/H1N1 virus revealed by steered molecular dynamics approach. J Chem Inf Model;50:2236–2247. 22. Mai B.K., Li M.S. (2011) Neuraminidase inhibitor R125489–a promising drug for treating influenza virus: steered molecular dynamics approach. Biochem Biophys Res Commun;410:688–691. 23. Lu H., Isralewitz B., Krammer A., Vogel V., Schulten K. (1998) Unfolding of titin immunoglobulin domains by steered molecular dynamics simulation. Biophys J;75:662–671. Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes
24. Colizzi F., Perozzo R., Scapozza L., Recanatini M., Cavalli A. (2010) Single-molecule pulling simulations can discern active from inactive enzyme inhibitors. J Am Chem Soc;132:7361–7371. 25. Nicolini P., Frezzato D., Gellini C., Bizzarri M., Chelli R. (2013) Toward quantitative estimates of binding affinities for protein-ligand systems involving large inhibitor compounds: a steered molecular dynamics simulation route. J Comput Chem;34:1561–1576. 26. Humphrey W., Dalke A., Schulten K. (1996) VMD: visual molecular dynamics. J Mol Graph;14:33–38, 27–8. 27. Duan Y., Wu C., Chowdhury S., Lee M.C., Xiong G., Zhang W., Yang R., Cieplak P., Luo R., Lee T., Caldwell J., Wang J., Kollman P. (2003) A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J Comput Chem;24:1999–2012. 28. Ryckaert J.P., Ciccotti G., Berendsen H.J.C. (1977) Numerical-integration of Cartesian equations of motion of a system with constraints – molecular-dynamics of N-alkanes. J Comput Phys;23:327–341. 29. Jiang W., Hardy D.J., Phillips J.C., Mackerell A.D. Jr, Schulten K., Roux B. (2011) High-performance scalable molecular dynamics simulations of a polarizable
Chem Biol Drug Des 2015
force field based on classical Drude oscillators in NAMD. J Phys Chem Lett;2:87–92. 30. Chen L.Y. (2015) Hybrid steered molecular dynamics approach to computing absolute binding free-energy of ligand-protein complexes—A brute force approach that is fast and accurate. J Chem Theory Comput;11:1928–1938. 31. Uhrinova S., Uhrin D., Powers H., Watt K., Zheleva D., Fischer P., McInnes C., Barlow P.N. (2005) Structure of free MDM2 N-terminal domain reveals conformational adjustments that accompany p53-binding. J Mol Biol;350:587–598. 32. Chen H.F., Luo R. (2007) Binding induced folding in p53-MDM2 complex. J Am Chem Soc;129:2930– 2937. 33. Chen J., Wang J., Zhu W., Li G. (2013) A computational analysis of binding modes and conformation changes of MDM2 induced by p53 and inhibitor bindings. J Comput Aided Mol Des;27:965–974.
Note a
Case DA, Darden TA, Cheatham TE III, Simmerling CL, Wang J, Duke RE, et al. (2008) AMBER10. San Francisco, CA: University of California, San Francisco.
9
Characterizing the Free-Energy Landscape of MDM2 Protein–Ligand Interactions by Steered Molecular Dynamics Simulations Guodong Hu1,2,*, Shicai Xu1,2 and Jihua Wang1,2,* 1
Shandong Provincial Key Laboratory of Functional Macromolecular Biophysics, Dezhou University, Dezhou, China 2 College of Physics and Electronic Information, Dezhou University, Dezhou, China *Corresponding authors: Guodong Hu, [email protected]; Jihua Wang, [email protected] Inhibition of p53–MDM2 interaction by small molecules is considered to be a promising approach to re-activate wild-type p53 for tumor suppression. Several inhibitors of the MDM2–p53 interaction were designed and studied by the experimental methods and the molecular dynamics simulation. However, the unbinding mechanism was still unclear. The steered molecular dynamics simulations combined with Brownian dynamics fluctuation–dissipation theorem were employed to obtain the free-energy landscape of unbinding between MDM2 and their four ligands. It was shown that compounds 4 and 8 dissociate faster than compounds 5 and 7. The absolute binding free energies for these four ligands are in close agreement with experimental results. The open movement of helix II and helix IV in the MDM2 protein-binding pocket upon unbinding is also consistent with experimental MDM2unbound conformation. We further found that different binding mechanisms among different ligands are associated with H-bond with Lys51 and Glu25. These mechanistic results may be useful for improving ligand design. Key words: MDM2 protein, steered molecular dynamics simulations, unbinding mechanism Received 12 April 2015, revised 20 May 2015 and accepted for publication 25 May 2015
The human version of mouse double minute protein 2 (MDM2) inhibits tumor suppressor activity of p53 through binding of its N-terminal domain to the transactivation domain of p53 (1–4). Disruption of the binding between MDM2 and p53 would enhance the ability of p53 to suppress tumor growth. As a result, MDM2 is an important anticancer target for small-molecule inhibitor design (5–8). ª 2015 John Wiley & Sons A/S. doi: 10.1111/cbdd.12598
Using a structure-based approach, Ding et al. (9) successfully designed several inhibitors of the MDM2–p53 interaction. The most potent ligand 8 has a Ki value of 3 nM binding to MDM2 and is 2000 times more potent than the natural p53 peptide (residues 13–19). These inhibitors were designed to mimic the amino acid residues in the short a-helical region of p53 (Phe19, Trp23, and Leu26) that interacts with MDM2. Molecular dynamics simulations coupled with free-energy calculations by molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) confirmed those ligands as residue mimicries although not all ligands mimicked the particular residue as designed, such as ligands 4 and 8 with a 2-morpholin-4yl-ethylamine group. However, the information about the unbinding of ligands from the binding pocket of MDM2 was still unclear (10). The purpose of this study was to examine the energy landscape of ligand–MDM2 interactions by steered molecular dynamics (SMD) simulations. SMD takes inspiration from single-molecule pulling experiments (11) and dissociates a complex structure by a pulling force (12,13). The non-equilibrium dynamics of the system under a pulling force can extract the free-energy difference between the bound state and the unbound state by measuring the work along the transition paths and map out the freeenergy landscape in terms of the potential of mean force (PMF) (14,15) with high precision and efficiency (16–18). SMD simulations have become widely used in studying many biochemical processes including the transportation of ions and organic molecules across membrane channels (16,19–22) and unfolding/folding mechanism of proteins (23,24) as well as the mechanisms of protein–ligand binding (12,17,25). In this study, we combined SMD with Brownian dynamics fluctuation–dissipation theorem (BD-FDT) to calculate the interaction energy landscape of MDM2 with four known MDM2 ligands (compounds 4, 5, 7, and 8, Figure 1) (9). The largest difference between ligand 8 and ligand 5 is the 2-morpholin-4yl-ethylamine group. The van der Waals energies, electrostatic energies, and hydrogen bonds were analyzed in detail along the pulling paths of SMD simulation. The binding affinities of these four ligands from energy landscape studies are in close agreement with experimental measurements. 1
Hu et al.
Figure 1: Molecular structures of four ligands which are labeled with their name as in reference (7) (compounds 4, 5, 7 and 8).
Model and Method Systems setup The starting structures for MDM2–ligand complexes without water molecules and counter ions were obtained from the last conformation of 3-nseconds equilibrated MD simulations of our previous work (10). To minimize steric
crashes during SMD simulations, all initial complex structures were rotated so that the orifice of MDM2’s binding pocket points to +z direction (Figure 2). The complexes were then solvated in a rectangular water box using VMD software (26). The minimum distances from the surface of complexes to the edges of boxes are 15 A in x, y, and –z directions and 25 A in the +z direction. Afterward, an
Figure 2: The orientations of ligands in the binding pocket of MDM2 from +z direction to the complex. The ligands are shown in stick and ball representations; the pulled atoms are shown in large balls.
2
Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes
appropriate number of counterions were added to neutralize the charge of systems. They were placed at the grids with the largest positive coulombic potentials around the complexes using the Tleap module implemented in a AMBER10 package . The standard AMBER force field for bioorganic systems (FF03) (27) was utilized to describe the proteins and water molecule (TIP3P). The parameters for ligands were obtained previously followed a standard AMBER procedure for force-field parameter and charge assignment through quantum calculations (10).
Molecular dynamics (MD) simulation and steered molecular dynamics (SMD) simulation The MD simulations and SMD simulations were performed using NAMD 2.8 package (15). We employed periodic boundary conditions in all three dimensions and utilized particle mesh Ewald (PME) method for treating long-range electrostatic interactions. The cutoff distance for longrange interactions was set to 12 A with a switching distance of 10 A. The time step was 2 fseconds for shortrange and 4 fseconds for long-range interactions. The covalent bonds of all hydrogen atoms were fixed to their equilibrium length by SHAKE algorithm (28). The temperature is maintained at 300 K using Langevin thermostat (29) with the damping constant chosen as 5/pseconds and pressure at 1 bar. To remove possible steric crashes between water molecules and complex structures, the energy of each system was minimized for 2000 steps, during which locations of protein and ligand atoms were restrained by a harmonic potential with a force constant of 0.2 kcal/(mol 9 A2). Afterward, 5-nseconds isothermal isobaric ensemble (NPT)-MD simulations were performed without any restraints. In our SMD, a ligand was pulled from a bound state to a dissociated state. The mass center of a ligand was initially set at the origin of co-ordinates. To maintain the location of MDM2 during the SMD simulations, we restrained the z-coordinates of several backbone carbon atoms of MDM2 that are far away from the binding pocket (Figure 3A). Four MDM2–ligand complex systems were equilibrated for 150 000 steps for each segment without pulling forces and for 50 000 steps with pulling forces.
Potential of mean force (PMF) from steered molecular dynamics (SMD) simulation Because the four ligands studied in this study have more than seventy atoms each, one center for a ligand is unlikely to sample all relevant events. As a result, each ligand is represented by two centers (Figure 2). Both centers were steered along z direction. The pulling speed was set at 0.01 A/pseconds at the z direction, and a ligand can move freely in x and y directions. To reduce the impact of pulling on MDM2, we divided the pulling path into 16 segChem Biol Drug Des 2015
A
B
Figure 3: (A) Pulling compound 8 from its bound state to an unbound state. MDM2 (a new cartoon and a surface representation), compound 8 (a ball representation), and the pulling path (the red line) are shown along with the displacements along z-axis and in xy plane as labeled. The restrained atoms are shown in a ball representation with gray color. (B) Potential of mean force for four ligands (compounds 4, 5, 7, and 8) as a function of the displacement from their binding site along the pulling path when two centers were steered away from the MDM2 protein.
ments along the z direction with 1 A for each segment. For each segment, we performed two types of SMD simulations: one for pulling away (denoted as forward) and one for pulling back to (denoted as reverse) the binding site. For each segment, two end states are denoted as A and B, respectively. Eight forward and reverse pulling paths between A and B states are sampled during which the work done to the system was recorded. The Gibbs freeenergy difference between state A and an intermediate state r (namely, the PMF) was given by the Brownian dynamics fluctuation–dissipation theorem (BD-FDT) (14,30) as follows: hexp½WA!r =2kB TiF DG ¼ kB T ln hexp½Wr!A =2kB TiR
where kB and T are the Boltzmann constant and the absolute temperature, respectively. The brackets with subscript F/R represent the statistical average over the forward/ reverse paths, respectively. WA?r is the work done to the 3
Hu et al.
Figure 4: Work done along the forward and reverse pulling paths. In the forward direction, the atoms of a ligand were pulled away from the binding site. In the reverse direction, the atoms of a ligand were pulled back to the binding site.
system along a forward path when the ligand is steered from A to r. Wr?A (=WB?A – WB?r) is the work done to the system along a reverse path when the ligand is steered from r to A. Data analysis We define that a ligand and MDM2 are in contact if an atom of MDM2 is within 7.0 A of any atoms of a ligand. The hydrogen bonds (H-bonds) are defined by the distance between acceptor and donor atoms ( 13 A, indicating that a 25- A water box for the +z direction is sufficient for the SMD simulations. Compounds 4 and 8, however, reached plateau at about 10 A, compared to 12 A for compounds 5 and 7; that is, the former can dissociate more quickly than the latter. Compound 5 reveals an interesting intermediate state at around 4–6 A. On the other hand, the PMF for compound 7 is made of two phases: a rapid increase phase from 0 to 3 A and a steady but slower increase phase after 3 A. Thus, the binding thermodynamics of four ligands are different. To understand different mechanisms for different compounds observed above, we calculated the interacting Table 1: The binding free energies for four complexesa Compound b
Figure 4 compares the work required for forward and reverse directions. The figure shows that the pulling process is largely irreversible on the basis of the following details. WA?z is not equal to Wz?A in most of the segments, especially in the first few segments of pulling. WA? z would be equal to Wz?A if the pulling was reversible. 4
MM-PBSA BD-FDT (this work) Expc
4
5
7
8
16.00 9.93 10.83
10.74 9.53 10.19
12.36 9.59 10.1
19.64 10.88 11.7
a
All values are give in kcal/mol. The binding free energies were calculated using the MM-PBSA method (8). c The experimental binding free energies were calculated using Ki (7). b
Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes A
B
C
Figure 5: The averaged energies as a function of ligand displacement along z-axis were (A) for van der Waals energies and (B) for electrostatic energies. (C) The averaged number of Hbonds between ligands and MDM2 as a function of ligand displacement along z-axis; the key values were labeled with ligand name and letters.
energy terms for the structures along the pulling path. For a given z value, van der Waals energies and the electrostatic energies were averaged over 16 structures extracted from eight forward SMD trajectories and eight reverse SMD trajectories. Figure 5A shows that the van der Waals energies increase smoothly as a function of the ligand displacement along the z-axis for four complexes. The magnitudes of van der Waals energies for four complexes at z = 0 have the same rank order as their experimental binding affinity values. Smooth increase of van der Waals energies as z increases indicates lack of steric hindrance for all ligands transiting from a abound state to an unbound state. However, large fluctuations of electrostatic energies are observed (Figure 5B). The minima for electrostatic Chem Biol Drug Des 2015
energies are located between 5 A to 8 A for compound 4, 8 A to 12 A for compound 5, and 8 A to 10 A for compound 7. To have a better understanding of large fluctuation in electrostatic energies, we calculated the number of H-bonds between a ligand and MDM2 (Figure 5C). Although the average number of protein–ligand H-bonds is mostly smaller than one, all the large dips in electrostatic energies for compounds 4, 5, and 7 correspond to a peak in the number of H-bonds. This indicates the importance of hydrogen bonding energies in overall electrostatic energies, despite that the number of H-bonds is mostly 0.5) with MDM2 until z > 8 A. The H-bonds were formed between its N3 H6 and the oxygen atoms of side chain of the flexible N-terminal residue Glu25 which allows a longer distance from compound 5 to MDM2 (Figure 7C). These H-bonds are semistable as their probability is about 0.5 (Figure 5C).
Unbinding mechanisms of compounds 7 and 8 Unlike compounds 4 and 5, the H-bonds between compounds 7/8 and MDM2 decay more slowly. The number of H-bond between MDM2 and compound 7 decreases smoothly from z = 0 to z = 6 A and then increases with a small peak at z = 8 A. For compound 8, its H-bond to MDM2 starts to decrease only at z = 6 A and becomes nearly zero at z > 8 A. Compounds 7 and 8 both had an H-bond to the oxygen atom of backbone of Leu54 (Figures 8B and 9) as with compounds 4 and 5. However, unlike compounds 4 and 5, geometry of compounds 7 and 8 allowed a large displacement as shown in Figure 9. Unlike compound 8, the atoms N3 H6 of ligand 7 can form a new H-bond with the oxygen atom of backbone of Lys51 at z = 8 A. Figure 9 shows the overlapped four structures extracted from the SMD simulation with the center of ligand at the position shown in Figure 5C. In summary, the H-bonds between nitrogen atom of a ligand and the backbone oxygen atom of Leu54 played 5
Hu et al.
B A
D
C
Figure 6: Schematic view of the key distances between ligand 4 and MDM2. MDM2 was shown in new cartoon representation and colored by secondary structure. Ligand was shown in stick representation and colored by name of atoms. The key atoms were shown in ball representation. The distance was labeled with dash line. The plots A, B, C, and D were structures for the center of ligand 4 at 0 A, 3.93 A, 5.36 A, and 7.11 A, respectively.
A
B
D
a key role for MDM2–ligand binding, consistent with previous work (9,10). Their stabilities are different for different compounds. After breaking the H-bond to Leu54, compounds 4 and 7 can form H-bond with Lys51, 6
C
Figure 7: Schematic view was similar to Figure 6 except for ligand 5. The plots A, B, C, and D were structures for the center of ligand 5 at 0 A, 3.06 A, 10.63 A, and 12.3 A, respectively.
whereas compound 5 forms H-bond with Glu25, but not compound 8. Thus, Lys51 and Glu25 are important for ligand dissociations from MDM2 for three of four ligands studied. Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes
Figure 8: Schematic view was similar to Figure 6 except for ligand 7. The plots A, B, C, and D were structures for the center of ligand 7 at 0 A, 3.74 A, 6.56 A, and 8.25 A, respectively.
A
B
C
D
Figure 9: Four structures of 8a (0 A), 8b (4.69 A), 8c (7.55 A), and 8d (8.54 A) were overlapped according to MDM2. The key section was enlarged and labeled.
The conformation of the binding cleft To study the effect of ligand dissociations on the dynamics and conformations in the binding pocket of the MDM2, RMSFs versus the residue number for MDM2 in the bound (z ≤ 1 A) and unbound states (15 A < z ≤ 16 A) were analyzed (Figure 10). The RMSFs of the most residues in MDM2 in the unbound state were larger than those in the bound state, the residues 41–56, and residues 91–100 that consist of loop regions and portions of helix II and helix IV. We calculated the averaged distances of four compounds between the center of helix II and helix IV in the bound and unbound states. The averaged fluctuation in the unbound state increases by about 0.43 A than the bound state, consistent with 0.5 A and 0.7 A increase upon p53 unbinding (31–33). Chem Biol Drug Des 2015
Conclusions Steered molecular dynamics combined with free-energy calculations is emerging as a promising tool in structuredbased drug design. In this study, SMD combined with BDFDT was used to study the interaction of four ligands with the MDM2 protein. The PMFs of these four ligands indicate that compounds 4 and 8 can reach an unbound state faster than compounds 5 and 7. This is likely because compounds 5 and 7 can have H-bonds with Lys51 at a large displacement from MDM2. The binding affinities from the PMF were in good agreement with the experimental data. This work provides the atomistic details of the unbinding mechanism, revealing transient interactions at accessory binding sites which can be exploited in-depth for inhibitor design for MDM2. 7
Hu et al.
Figure 10: Root mean square fluctuation of backbone atoms versus residue number of MDM2. The regions of four helices of MDM2 were labeled.
Acknowledgments This work was partially supported by funding from the National Natural Science Foundation of China (61271378 and 11447004), the Natural Science Foundation of Shandong Province (ZR2014JL006 and ZR2012CL09), and A Project of Shandong Province Higher Educational Science and Technology Program (J14LJ05). We would like to thank Yaoqi Zhou for his scientific support.
Conflict of Interests The authors declare that they have no conflict of interest.
References 1. Vogelstein B., Lane D., Levine A.J. (2000) Surfing the p53 network. Nature;408:307–310. 2. Levine A.J. (1997) p53, the cellular gatekeeper for growth and division. Cell;88:323–331. 3. Zhang X.P., Liu F., Cheng Z., Wang W. (2009) Cell fate decision mediated by p53 pulses. Proc Natl Acad Sci U S A;106:12245–12250. 4. Zhang X.P., Liu F., Wang W. (2011) Two-phase dynamics of p53 in the DNA damage response. Proc Natl Acad Sci U S A;108:8990–8995. 5. Vassilev L.T., Vu B.T., Graves B., Carvajal D., Podlaski F., Filipovic Z., Kong N., Kammlott U., Lukacs C., Klein C. (2004) In vivo activation of the p53 pathway by small-molecule antagonists of MDM2. Science;303:844–848. 6. Ding K., Lu Y., Nikolovska-Coleska Z., Qiu S., Ding Y., Gao W., Stuckey J., Krajewski K., Roller P.P., Tomita Y., Parrish D.A., Deschamps J.R., Wang S. (2005) Structure-based design of potent non-peptide MDM2 inhibitors. J Am Chem Soc;127:10130–10131. 7. Vassilev L.T. (2005) p53 activation by small molecules: application in oncology. J Med Chem;48:4491–4499. 8
8. Hang Y., Gui-in L., Hyung Soon P., Gregory A.P., Johanna M.R., Said M.S., Andrew D.H. (2005) Terphenyl-based helical mimetics that disrupt the p53/HDM2 interaction13. Angew Chem Int Ed;44:2704–2707. 9. Ding K., Lu Y., Nikolovska-Coleska Z., Wang G., Qiu S., Shangary S., Gao W., Qin D., Stuckey J., Krajewski K., Roller P.P., Wang S. (2006) Structure-based design of spiro-oxindoles as potent, specific small-molecule inhibitors of the MDM2-p53 interaction. J Med Chem;49:3432–3435. 10. Hu G., Wang D., Liu X., Zhang Q. (2010) A computational analysis of the binding model of MDM2 with inhibitors. J Comput Aided Mol Des;24:687–697. 11. Grubmuller H., Heymann B., Tavan P. (1996) Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force. Science;271:997–999. 12. Patel J.S., Berteotti A., Ronsisvalle S., Rocchia W., Cavalli A. (2014) Steered molecular dynamics simulations for studying protein-ligand interaction in cyclindependent kinase 5. J Chem Inf Model;54:470–480. 13. Li M.S., Mai B.K. (2012) Steered molecular dynamics – a promising tool for drug design. Curr Bioinform;7:342. 14. Chen L.Y. (2008) Nonequilibrium fluctuation-dissipation theorem of Brownian dynamics. J Chem Phys;129:144113. 15. Allen T.W., Andersen O.S., Roux B. (2006) Molecular dynamics – potential of mean force calculations as a tool for understanding ion permeation and selectivity in narrow channels. Biophys Chem;124:251–267. 16. Chen L.Y. (2013) Glycerol modulates water permeation through Escherichia coli aquaglyceroporin GlpF. Biochim Biophys Acta;1828:1786–1793. 17. Chen L.Y. (2011) Exploring the free-energy landscapes of biological systems with steered molecular dynamics. Phys Chem Chem Phys;13:6176–6183. 18. Hu G., Chen L.Y. (2010) In silico experiments of single-chain antibody fragment against drugs of abuse. Biophys Chem;153:97–103. 19. Giorgino T., De Fabritiis G. (2011) A high-throughput steered molecular dynamics study on the free energy profile of ion permeation through gramicidin A. J Chem Theory Comput;7:1943. 20. Henin J., Tajkhorshid E., Schulten K., Chipot C. (2008) Diffusion of glycerol through Escherichia coli aquaglyceroporin GlpF. Biophys J;94:832–839. 21. Mai B.K., Viet M.H., Li M.S. (2010) Top leads for swine influenza A/H1N1 virus revealed by steered molecular dynamics approach. J Chem Inf Model;50:2236–2247. 22. Mai B.K., Li M.S. (2011) Neuraminidase inhibitor R125489–a promising drug for treating influenza virus: steered molecular dynamics approach. Biochem Biophys Res Commun;410:688–691. 23. Lu H., Isralewitz B., Krammer A., Vogel V., Schulten K. (1998) Unfolding of titin immunoglobulin domains by steered molecular dynamics simulation. Biophys J;75:662–671. Chem Biol Drug Des 2015
Unbinding of MDM2-Ligand Complexes
24. Colizzi F., Perozzo R., Scapozza L., Recanatini M., Cavalli A. (2010) Single-molecule pulling simulations can discern active from inactive enzyme inhibitors. J Am Chem Soc;132:7361–7371. 25. Nicolini P., Frezzato D., Gellini C., Bizzarri M., Chelli R. (2013) Toward quantitative estimates of binding affinities for protein-ligand systems involving large inhibitor compounds: a steered molecular dynamics simulation route. J Comput Chem;34:1561–1576. 26. Humphrey W., Dalke A., Schulten K. (1996) VMD: visual molecular dynamics. J Mol Graph;14:33–38, 27–8. 27. Duan Y., Wu C., Chowdhury S., Lee M.C., Xiong G., Zhang W., Yang R., Cieplak P., Luo R., Lee T., Caldwell J., Wang J., Kollman P. (2003) A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J Comput Chem;24:1999–2012. 28. Ryckaert J.P., Ciccotti G., Berendsen H.J.C. (1977) Numerical-integration of Cartesian equations of motion of a system with constraints – molecular-dynamics of N-alkanes. J Comput Phys;23:327–341. 29. Jiang W., Hardy D.J., Phillips J.C., Mackerell A.D. Jr, Schulten K., Roux B. (2011) High-performance scalable molecular dynamics simulations of a polarizable
Chem Biol Drug Des 2015
force field based on classical Drude oscillators in NAMD. J Phys Chem Lett;2:87–92. 30. Chen L.Y. (2015) Hybrid steered molecular dynamics approach to computing absolute binding free-energy of ligand-protein complexes—A brute force approach that is fast and accurate. J Chem Theory Comput;11:1928–1938. 31. Uhrinova S., Uhrin D., Powers H., Watt K., Zheleva D., Fischer P., McInnes C., Barlow P.N. (2005) Structure of free MDM2 N-terminal domain reveals conformational adjustments that accompany p53-binding. J Mol Biol;350:587–598. 32. Chen H.F., Luo R. (2007) Binding induced folding in p53-MDM2 complex. J Am Chem Soc;129:2930– 2937. 33. Chen J., Wang J., Zhu W., Li G. (2013) A computational analysis of binding modes and conformation changes of MDM2 induced by p53 and inhibitor bindings. J Comput Aided Mol Des;27:965–974.
Note a
Case DA, Darden TA, Cheatham TE III, Simmerling CL, Wang J, Duke RE, et al. (2008) AMBER10. San Francisco, CA: University of California, San Francisco.
9
