Determination of the Structure of Human Phosphodiesterase-2 in a

Feb 9, 2009 - Fundamental Reaction Mechanism and Free Energy Profile for (−)-Cocaine Hydrolysis Catalyzed by Cocaine Esterase. Junjun Liu , Adel ...
2 downloads 0 Views 3MB Size
Download PDF
2896

J. Phys. Chem. B 2009, 113, 2896–2908

Determination of the Structure of Human Phosphodiesterase-2 in a Bound State and Its Binding with Inhibitors by Molecular Modeling, Docking, and Dynamics Simulation Adel Hamza and Chang-Guo Zhan* Department of Pharmaceutical Sciences, College of Pharmacy, UniVersity of Kentucky, 725 Rose Street, Lexington, Kentucky 40536 ReceiVed: September 16, 2008; ReVised Manuscript ReceiVed: NoVember 13, 2008

Although an X-ray crystal structure of the catalytic domain of human cyclic nucleotide phosphodiesterase-2 (PDE2A) was reported in the literature, the shape of the binding pocket is not suitable for binding with its known potent inhibitors. Extensive molecular modeling, docking, and dynamics simulations in the present study have demonstrated for the first time that the structure of PDE2A catalytic domain may exist in two different forms corresponding to the inhibitor-bound and unbound states of the enzyme. The structural change from the unbound state to the bound state leads to a substantial variation in the size of the pocket but does not affect the general structural feature of the catalytic site. The flexible binding pocket and conserved structural feature of the catalytic site lead us to better understand why this enzyme can catalyze hydrolysis of two different intercellular second messengers. It has been demonstrated that the available X-ray crystal structure of PDE2A was in the unbound state, explaining why it is not suitable for molecular docking studies on the enzyme-inhibitor binding. We have developed a reasonable 3D model of the PDE2A structure in the bound state and determined the detailed binding modes and binding free energies for PDE2A binding with its known potent inhibitors. The calculated binding free energies are in good agreement with available experimental data. The general structural insights, PDE2A model in the bound state, and detailed PDE2A-inhibitor binding structures obtained in this study will be valuable for future rational design of novel, potent inhibitors of PDE2A as therapeutic agents. Introduction The cyclic adenosine 3′,5′-monophosphate (cAMP) or cyclic guanosine 3′,5′-monophosphate (cGMP) are involved in a number of intracellular processes, such as signal transduction, gene transcription, activation of kinases, and regulation of channel function, as the second messengers.1 High intracellular level of intracellular second messenger cAMP or cGMP can be achieved either by activation of the adenylyl/guanylyl cyclases or by inhibition of the cyclic nucleotide-degrading enzymes phosphodiesterases (PDEs). This superfamily of distinct enzymes metabolizes cAMP and cGMP to their biologically inactive nucleotide 5′-monophosphates AMP and GMP. PDEs are ubiquitous enzymes that control the intracellular concentrations of the cyclic nucleotides by catalyzing their hydrolysis,2,3 and these enzymes are essential for cyclic nucleotide signaling since they represent, at least for eukaryotic cells, essentially the only way for a cell to terminate a cyclic nucleotide signal. In addition, PDEs are largely responsible for confining a cyclic nucleotide signal to particular locations and to prevent its diffusion throughout the cell.4-6 The human genome codes for 11 PDE families resulting in a PDE proteome of about 60 PDE isoenzymes.2,7 The major families involved in cGMP hydrolysis are PDE1, PDE2, PDE5, PDE6, and PDE9.8,9 PDE2 is cGMP-stimulated and metabolizes both cGMP and cAMP, although its affinity for cGMP is slightly higher than that for cAMP.10 High PDE2 activity can be found in the heart11 and brain.10 Lower expression of PDE2A3 mRNA was found in the lung, placenta, liver, skeletal muscle, kidney, and pancreas.10,12 PDE2 exists in a homodimer of two 105 kDa * To whom correspondence should be addressed. Tel: 859-323-3943. Fax: 859-323-3575. E-mail: [email protected].

subunits. Each subunit contains 2 GAF domains (found in cGMP-selective phosphodiesterases, adenylyl cyclases, and formate hydrogen lyase transcriptional activator (FhlA)), which allows either the dimerization or binding of cGMP.13 In contrast to other PDE isoenzyme families, PDE2 is encoded by one single gene and only one human splice variant has been described so far. The catalytic domains of all PDE family members share considerable homology in their amino acid sequences, and their three-dimensional (3D) structures appear to be very similar. Despite the overall similarity of their catalytic domains, the individual PDEs exhibit characteristic selectivity for cAMP and/ or cGMP14 as their substrates, and each displays a unique inhibitor specificity profile. Besides their interest for basic cell biology and biochemistry, the human PDEs have become interesting targets for drug discovery and development. Familyand tissue-specific PDE inhibitors constitute a growing class of pharmaceutical compounds that find applications for a broad spectrum of maladies.15-17 Thus, design of molecules that selectively inhibit specific PDEs is of considerable biological significance and will likely lead to novel PDE inhibitors for important clinical applications. Within the past few years, X-ray crystal structures for the catalytic domains from representative members of different PDE families have been solved. Several reviews18-20 have appeared describing the details of the structure-function relationships discovered from these studies. In addition, the structure of the catalytic domain and of the regulatory tandem GAF domains of PDE2A has also been reported.13,21 However, no highresolution X-ray crystal structure for any PDE holoenzyme has

10.1021/jp8082612 CCC: $40.75  2009 American Chemical Society Published on Web 02/09/2009

Structure of PDE2 in a Bound State been reported, so we know little about the molecular details of how the regulatory domains influence catalysis. The overall folds and functional/structural elements of the catalytic domains of different PDEs are very similar despite that there is only about 25% to 35% amino acid sequence identity between the catalytic domains of different PDE proteins. Thus, these catalytic domains contain three subdomains composed largely of 16 helices. The active site is formed at the junction of the helices by residues that are highly conserved among all the PDEs. At the bottom of the substrate-binding pocket are two divalent metal binding sites. The two metal (zinc and magnesium) ions are coordinated by residues located on each of the three different domains. The metal-binding site that binds zinc has two histidine residues and two aspartic acid residues that are conserved among all PDEs studied to date. These residues form part of the signature recognition sequence for cyclic nucleotide PDEs. Moreover, several of reported X-ray crystal structures of PDEs have inhibitors bound inside the binding pocket. It appears from these structural studies that, although PDE inhibitors bind at the active site, three different binding modes have been found to date.18,20 Indeed, X-ray crystal structures for PDE4B, PDE4D, and PDE5A binding with inhibitors revealed that the compounds interact with the enzymes either through hydrogen bonds with residues involved in the nucleotide binding, through interactions with hydrophobic residues lining the active site channel, or with metal ions mediated through water.18 The variety of mechanisms of inhibitor binding has encouraged medicinal chemists in their quest for increasingly selective drugs. Elucidation of more PDE structures in the absence or presence of inhibitors should facilitate the design and discovery of more potent and selective PDE inhibitors in the future. The recently reported X-ray crystal structure of the human PDE2A catalytic domain was solved in an apo form or unbound state.21 The phosphate group found in the binding pocket is quite small to be considered as an inhibitor. The authors21 also proposed the possible binding mode of the known ligands Rolipram (a PDE4 inhibitor), erythro-9-(2-hydroxy-3-nonyl) adenine (EHNA), and AMP in the PDE2A binding site, without carrying out any computational modeling (e.g., molecular docking, energy minimization, or simulation). EHNA is the only PDE2-specific inhibitor that is commercially available. EHNA has been demonstrated to specifically act on PDE2 by inhibiting cGMP activation of PDE2 with an IC50 value of ∼1 µM and an at least 50-fold selectivity over other PDEs.22-24 The core structure of EHNA resembles cAMP but differentiates in the fact that EHNA has a bulky hydrophobic carbon side chain replacing the phosphoribose moiety in cAMP. More potent PDE2 inhibitors have been reported recently. The most promising potent PDE2 inhibitors known so far include ND7001, a benzodiazepinone derivate from Neuro3D Co.,25 with an IC50 value of ∼50 nM and compound Bay607550 (Figure 1) with an IC50 value of ∼4.7 nM for human PDE2A.23 These potent inhibitors of PDE2 can be used as the lead for further design of new compounds that are more potent for human PDE2. For rational design of more potent inhibitors of PDE2 starting from these available lead compounds, we first need to understand how these potent inhibitors bind with PDE2. However, it is unclear how these important ligands bind with PDE2 in the active site and, to the best of our knowledge, no computational simulation on PDE2 inhibitor has been reported in the literature so far. A problem for computational studies of PDE2-inhibitor binding is that the active site of X-ray crystal structure available for PDE2 was in a collapsed form (unbound

J. Phys. Chem. B, Vol. 113, No. 9, 2009 2897

Figure 1. 2D structures of Bay60-7550 (left) and benzodiazepinone derivative (ND7001) (right) compounds.

state) of the enzyme without a meaningful ligand. Molecular docking studies cannot be realized by simply using the available X-ray crystal structure. In the present study, we aimed to understand how human PDE2A binds with its potent inhibitors and, thus, to establish a solid structural base for future rational design of more potent PDE2A inhibitors as potential therapeutic agents. For this purpose, we first developed a reasonable 3D model of PDE2A representing the enzyme structure in a bound state through computational modeling using the high-resolution X-ray crystal structure of PDE5 as a template. Using the developed 3D structural model, we were able to carry out molecular docking, molecular dynamics (MD) simulations, and binding free energy calculations to explore possible binding modes of the potent PDE2A inhibitors and to identify key amino acids involved in the PDE2-inhibitor binding. All of the computational results reveal valuable insights into the microscopic binding of PDE2 with inhibitors and the structure-activity correlation. Based on the determined binding structures, the calculated binding free energies are in good agreement with available experimental data, suggesting that the computationally determined binding modes are reasonable. The structural and energetic insights obtained in this study are expected to be valuable for rational design of novel, more potent inhibitors of human PDE2. Methods Construction of Initial 3D Models of PDE2A. The amino acid sequence of human phosphodiesterase 2A (PDE2A) [AC: O00408] was generated from the GenBank database.26 The search for sequence similarity with several members of the PDE superfamily available in the Protein Data Bank (PDB)27 was performed by using BLAST program.28 The crystal structure of the human phosphodiesterase 5A (PDE5A) [ID 1XOZ, 1.37 Å resolution]18 was used as the template to build the initial PDE2A model. The multiple-sequence alignment was performed using Homology module of InsightII software (Accelrys, Inc.). The 3D models of the PDE2A (denoted by PDE2AHomology) were generated using the automated homology modeling tool Modeller (version 9v1)29-31 with the default parameters, while the remaining loops were built de novo and refined to access proper folding with minimum steric clash by using the loop subroutine of Modeller. Atomic coordinates of the side chains were added from the standard residues library of Modeller. The Modeller is a well-known comparative modeling tool, which generates a refined 3D homology model of a protein sequence automatically and rapidly, based on a given sequence alignment to a known 3D protein structure template. It employs probability density functions (PDFs) as the spatial restraints rather than energy.29-31 The Modeller generates a large number of templatederived restraints to force the model toward the structure of the template and converges to the best possible structure by

2898 J. Phys. Chem. B, Vol. 113, No. 9, 2009 simultaneously satisfying a network of spatial restraints and molecular geometry using the CHARMM force field.32 The optimization process to generate the 3D model consists of applying the variable target function as well as conjugate gradients and molecular dynamics (MD) with simulated annealing. Thus, to construct a preliminary 3D structural model as an input of further refinement, a set of 25 conformations of the 3D structure were generated using the fast simulatedannealing procedure implemented in the Modeller. It is wellknown that MD simulations are valuable in the improvement of medium resolution (∼4 Å) structures and homology models. As well as restricted MD has been used in some solved structure refinement programs,33 it has also been used in the completion/ refinement of homology models.34 Validation of the Models. In general, every homology model contains errors in the initial structure of the model. The number of errors for a given method is mainly dependent on the percentage sequence identity between the template and target and on the number of errors in the template itself. An essential step in homology modeling process is to verify/validate the model. We used several steps to estimate possible errors in our 3D models. First, the top 10 models obtained with the corresponding lowest probability density functions (PDFs) were selected. The stereochemical quality of these structures was evaluated to determine whether the bond lengths and bond angles are within their normal ranges, and whether there are lots of bumps in the models (corresponding to a high van der Waals energy). This was done by calculating G-factor using Procheck program.35,36 Second, the atomic contact analysis was also performed by using the program to identify bad packing of side chain atoms or unusual residue contacts. Third, ProsaII program37 was used to check the native folding of the modeled structures by calculating the z-score, which is the score of pair and surface energies based on mean force potentials of the model. In order to assess the reliability of each model, the corresponding energy graphs were compared with the structure of PDE5 template and the resolved PDE2A crystal structure. Then, the normality indices (z-score) and the atomic contact analysis that describe how well a given characteristic of the model resembles the same characteristic in the template structure were determined by using Whatif program.38 Finally, the 3D model with the best score (in terms of G-factor and z-score) was considered to the best initial structure for further refinement as described below. Model Refinement. Finally, the built initial PDE2A model was solvated in water and subjected to a short run of energy minimization (3000 cycles) to relieve possibly unfavorable steric interactions and to optimize the stereochemistry (see below for MD simulation procedure). This refined 3D model of PDE2A was further evaluated periodically for its stereochemical quality as well as its residue packing and atomic contacts as described above. Molecular Docking. Molecular geometries of Bay60-7550 and ND7001 were optimized by performing ab initio electronic structure calculations using the Gaussian03 program39 at the HF/ 6-31G* level. The optimized geometries were used to calculate the electrostatic potentials on the molecular surfaces at the same HF/6-31G* level. The calculated electrostatic potentials were used to determine the partial atomic charges by using the standard restrained electrostatic potential (RESP) fitting procedure.40 The determined RESP charges were used for the calculation of electrostatic energies in the docking and MD simulation processes. The missing force field parameters for the ligands were taken from the General Amber Force Field

Hamza and Zhan (GAFF) implemented in Amber9 program.41 The molecular docking for each ligand binding was carried out in the same way as we recently did for studying other protein-ligand binding systems.42-44 Docking calculations were performed on the ligands in the PDE2A binding site using the “automatic docking” Affinity module of the InsightII package (Accelrys, Inc.). The Affinity methodology uses a combination of Monte Carlo type and simulated Annealing procedures to dock the guest molecule (ligand) to the host (receptor).45 A key feature is that the “bulk” of the receptor, defined as atoms not in the binding (active) site specified, is held rigid during the docking process, while the binding site atoms and ligand atoms are allowed to move. During the initial docking calculation process, a roughly docked PDE2A-ligand complex was prepared. Thus, the first step of building the PDE2A-ligand complex was to dock the ligand to PDE2A by virtue of their geometric complementarity. We aimed to find where the ligand could be inserted most comfortably. The docking procedure consists of superimposing the PDE2AHomology (3D model) with the PDE2AX-ray-Rolipram complex proposed by Iffland et al.,21 and then the scaffold of the ligand was initially positioned and superimposed to the scaffold of the Rolipram to maximize the overlap between the two core structures. Thereafter, the other parts of the ligand were positioned in an approximate orientation on the PDE2A active site to make appropriate interactions without steric clashes between the ligand and the residues side chain of the pocket. Thus, we obtained a reasonable starting structure for further modeling studies. The Affinity/InsightII module was used to energy-minimize the starting structure. This process was designed to remove bad contacts and poor internal geometry in the initial structure and to obtain a reasonable starting point for subsequent search. Finally, the binding site was defined as a sphere centered on the initial position of the ligand with an approximately 10 Å radius around the catalytic active site of PDE2A. The internal volume of the PDE2A active site and molecular volumes of ligands were calculated with the MOLCAD module of SYBYL 7.0 (Tripos Inc., USA) using the Fast Connolly Channel algorithm with a probe size of 1.4 Å. The ligand was then moved by a random combination of translational, rotational, and torsional changes. The random movement of the ligand represents both the conformational space of the ligand and its orientation with respect to PDE2A. This procedure has the advantage of overcoming any energy barrier on the potential energy surface. For each resulting randomly moved structures, the energy was evaluated and compared to that of the previously energy-minimized structure. If the calculated energy change was within the specified (default) energy tolerance parameter, it was considered to have passed this first step and the structure was energy minimized. The second step was for fine tuning the docking. Whether the finally energy-minimized structure was accepted or rejected was based on the Metropolis energy criterion as implemented in the software and its similarity to structures found before. The Metropolis criterion was found to be best suited for finding a very small number of docked structures with very low energies (∼100 lowest-energy minimized structures were kept). We have further refined these PDE2A-ligand binding structures with the Simulated Annealing procedure (starting temperature at 360 K) and continued for 150 ps MD simulations at T ) 298.15 K. The energy minimization was performed by using the steepest descent algorithm first until the maximum energy derivative was smaller than 4 kcal mol-1 Å-1 and then

Structure of PDE2 in a Bound State

J. Phys. Chem. B, Vol. 113, No. 9, 2009 2899

Figure 2. (Left) Ribbon view of the clustered ND7001 binding poses (two major clusters) in the active site of the X-ray crystal structure of PDE2A. (Right) View of the most important residues close to the major clusters. The ligand conformations in the same cluster are displayed with the same color. Hydrogen bond is highlighted in dashed yellow line.

using the conjugated gradient algorithm until the maximum energy derivative was smaller than 0.001 kcal mol-1 Å-1. During the energy minimization and MD simulation, only the ligand and residues side chain in the binding pocket were kept free to move. The conformational searches (docking poses), energy minimization, and MD simulation for these processes were performed by using the Amber force field implemented in the Discover_3/InsightII calculation engine. The nonbonded interaction cutoff and the dielectric constant were set up to group based (20 Å cutoff distance) and distance dependent (ε ) 4r)46-48 to mimic the solvent environment, respectively. The MD simulation was performed with a time step of 1 fs. Cluster Analysis. The docked structures were examined, and the most favored cluster in the distribution of ligand poses was identified using a method similar to the one described by Gatchell et al.49,50 Briefly, for every docked pose, the number of ligand neighbors within a threshold root mean squared deviation (rmsd) was determined. The pose with the maximum number of neighbors, or the most populated cluster, was identified as the most favored pose or conformation. For the studied compounds, each docked pose was clustered by 0.8 Å of rmsd criterions. Single conformations or poses disconnected from the rest of the ensemble were considered as unlikely conformations and hence were identified and removed. After the initial criterion was satisfied, the second step was an examination of the different interactions that could be formed between the protein and the core structure of the ligand. In this way, the best cluster was selected based on the following criteria: lowest interaction energy (sum of the electrostatic and van der Waals interaction energy terms) of the PDE2A-ligand complex, hydrogen bonds formed between Gln859 side chain and the amide group of the ligand, and/or interaction with Phe862 side chain. Finally, the ligand conformation from the best cluster with the lowest interaction energy was subjected to energy minimizations and MD simulations on the explicitly solvated protein-ligand complexes (see below for the MD simulations) that can better account for the solvent effects on the protein-ligand binding. Molecular Dynamics Simulation Procedure. The initial geometry of PDE2A-ligand structure was neutralized by adding appropriate sodium counterions and was solvated in a rectangular box of TIP3P water molecules51 with a minimum solute-wall distance of 10 Å. Prior to MD simulations, the

protein (or protein-ligand complex) was kept fixed with a constraint of 500 kcal mol-1 Å-2, and we just energy-minimized the positions of the water molecules. Then, the entire solvated system was fully energy-minimized without any constraint. The MD simulation was performed by using the Sander module of Amber9 program41 in a way similar to what we did for other protein-ligand systems.52 The MD simulation was run by using the energy-minimized structure as the starting structure, and the particle mesh Ewald (PME) algorithm was used for dealing with long-range interactions.53 Each of the solvated systems was carefully equilibrated before a sufficiently long MD simulation in room temperature. The MD simulations were performed with a periodic boundary condition in the NPT ensemble at T ) 298.15 K with Berendsen temperature coupling54 and constant pressure (P ) 1 atm) with the isotropic molecule-based scaling. The time step of the simulation was 2.0 fs with a cutoff of 12 Å for the nonbonded interactions, and SHAKE algorithm was employed to keep all bonds involving hydrogen atoms rigid.55 Constant volume was carried out for 50 ps, during which the temperature was raised from 10 to 300 K. Then, a long production run of constant-pressure MD was carried out at 300 K. The first ∼2 ns of MD simulations on the PDE2A-ligand complexes were performed without any constraint, whereas a harmonic constraint of 50 kcal mol-1 Å-2 was applied on the R-carbon atoms during the last ∼1.6 ns MD simulations. We used the constraint for a couple of reasons. First, the MD simulation on the free PDE2A structure (see below for the results) revealed that the positions of the R-carbon atoms in the simulated free protein were very stable. In addition, it has been demonstrated42,43,56 that using certain constraint on R-carbon atoms of proteins could help to decrease the unnecessary fluctuations of the binding free energies calculated using different snapshots (see below for the binding free energy calculations). The atomic coordinates of the simulated complex were saved every 1 ps. The final structure of the complex was obtained as the average of the last 500 ps of MD minimized with the CG method until a convergence of 0.05 kcal mol-1 Å-1. The obtained stable MD trajectory was used to estimate the binding free energy (∆Gbind) by using the molecular mechanics/ Poisson-Boltzmann surface area (MM-PBSA) free energy calculation method.57 200 snapshots were used to perform the

2900 J. Phys. Chem. B, Vol. 113, No. 9, 2009 MM-PBSA calculation. Our MM-PBSA calculation for each snapshot was carried out in the same way as we did for other protein-ligand systems.42 The finally calculated binding free energy was taken as the average of the ∆Gbind values with the 200 snapshots. We note that the MM-PBSA binding energy method ignores the deformation energies of the protein and ligand themselves after the binding, which is based on an implicit assumption/approximation that the overall deformation energies are similar for all of the protein-ligand binding complexes and may be canceled out with other systematic approximations. In fact, the MM-PBSA method has been used to calculate the binding free energies of protein-protein, protein-nucleic acid, and protein-inhibitor interactions.57 A number of recent reports have demonstrated that this method is successful in predicting the protein-inhibitor binding free energies.58,59 Molecular Modeling of PDE2AX-rayfbound state. We also constructed another 3D model of the inhibitor-bound state of PDE2A starting from the X-ray crystal structure of PDE2A (the unbound state). The bound state structure modeled in this way is denoted by PDE2AX-rayfbound state. The PDE2AX-rayfbound state structure was obtained by removing the C-terminal fragment (residues 902 to 915) and rebuilding the residues 838 to 855 (loop LR14R15) of the resolved X-ray structure of PDE2A. The Modeler and Procheck software packages were also used to build and assess the stereochemical quality of the PDE2AX-rayfbound state structure. Then, the initial conformation of the loop was optimized and checked for its stability by submitting the PDE2AX-rayfbound state model structure to 1.7 ns of MD simulation in water system. The full protein was first relaxed during 500 ps of an unconstrained MD simulation then during the following 1.2 ns, only the residues of the LR14R15 loop and the side chain of Tyr827, Phe831, Ile855, and Leu858 residues in the upper binding pocket (residues close to loop) were free to move, while the rest of the protein was kept frozen (see MD simulation procedure). The average of the PDE2AX-rayfbound state structure was obtained by superimposing the snapshot structures collected every 200 ps of the MD simulation. Results and Discussion Binding Modes of ND7001 and Bay60-7550 with the X-ray Structure of PDE2A. Predicting reasonable binding modes of ligands in the protein binding pocket is a critical step when molecular docking simulation is used for structure-based drug design. At present, it is well-known that the binding of substrate (cGMP), or known inhibitors containing an pyrimidinone group, with PDE enzymes18,20 involves the formation of a bidentate hydrogen bond between the amide moiety of the pyrazolopyrimidinone group and the γ-amide group of Gln859 residue. This residue is completely conserved in all the PDE families, and an equivalent hydrogen bond was observed in previously reported PDE4/5-ligand complex structures.60-62 In addition to this interaction, the ligand is also stabilized by a π-π stacking interaction with the phenyl group of Phe862 residue. These interactions guided the ligand toward its binding site, leading to interactions with key residues in the PDE binding site. Moreover, these interactions are very important for understanding the catalytic behavior of the enzymes against substrate cGMP or cAMP.18,20 Thus, the refined 100 docked poses of ND7001 that passed the criteria described previously were clustered into two major and several minor binding modes. Snapshots for the two major clusters are represented in Figure 2. The centers of the two most populated ND7001 clusters are in the same location but with

Hamza and Zhan different orientations in the PDE2AX-ray binding pocket. In both clusters, the scaffold of the ligand is parallel to the phenyl ring of Phe862 residue and close to the amide group of Gln859 side chain. Although the orientation of the ligand in one cluster shows a monodentate hydrogen bond between the carbonyl amide group of the ligand and Gln859 side chain, none of the molecular conformations in the two clusters are stabilized by a π-π interaction with Phe862 side chain as observed in other PDE-ligand X-ray complexes.60-62 It seems like that the shape and the volume of the binding cavity of PDE2A cannot allow a reasonable pose of binding for ND7001 to be in a desirable deeper position. Moreover, manual superimposition of the Rolipram and EHNA structures onto ND7001 shows a net difference in size. Indeed, the molecular volumes are estimated to be ∼1076 Å3 for ND7001, ∼850 Å3 for Rolipram, and ∼880 Å3 for EHNA. Furthermore, an examination of the PDE2A binding pocket with the clustered ligand, highlighted several obvious steric obstructions that could impede the ligand binding. Indeed, during the docking and MD simulation, the residues in the upper domain (consisting of Ile826, Tyr827, Phe830, Phe831, Met847, Ala853, Ile855, and Leu858), more specifically Tyr827 and Phe830 side chains, appear to prevent a favorable binding of the ligand to Gln859 residue. To assess further the reliability of the docked ND7001 structures, the conformation with the lowest interaction energy from the two major clusters was refined by MD simulation in water. Thus, the MM-PBSA-calculated binding energies (without entropic contributions) for the refined ND7001 binding structures were found to be about -11 and -6 kcal/mol, respectively; including entropic contributions would make the calculated binding free energy values much higher and even become positive values. According to the high binding energy values calculated, the predicted conformation of ND7001 in the binding pocket clearly indicates that ND7001 does not fit the PDE2AX-ray active site pocket. We were unable to find a reasonable binding pose including the favorable interactions with Gln859 and Phe862 residues by directly using the X-ray crystal structure PDE2AX-ray for molecular docking studies. It should be pointed out that being unable to find a reasonable binding pose within nanoseconds of the MD simulation does not necessarily mean that a reasonable binding mode between the protein and ligand does not exist. This is because the large backbone change of a large protein like PDE2 is expected to be slow and will need a much longer time. In principle, an MD simulation could eventually lead to a reasonable protein-ligand binding structure, if the MD simulation could practically run for a sufficiently long time (e.g., milliseconds). We can only conclude from the MD simulation for nanoseconds that the backbone structure of PDE2AX-ray is not suitable for binding with the ligands. The similar molecular docking was also performed to examine PDE2AX-ray binding with Bay60-7550. First, the pyrazolopyrimidinone core of the ligand was positioned facing the phenyl ring of Phe862 residue and front of Gln859 side chain. This approximate pose allowed the Bay60-7550 scaffold to make favorable interaction with these two essential residues by establishing hydrogen bonds with amide group of Gln859 and a hydrophobic interaction with Phe862 side chain. The second step was to move and orient the two alkyl phenyl moieties into the active site pocket. After several rounds of docking tests, we finally concluded that it would not be possible to make an appropriate pose for having the alkyl metoxyphenyl chain in the pocket, due essentially to the presence of steric clashes with

Structure of PDE2 in a Bound State

J. Phys. Chem. B, Vol. 113, No. 9, 2009 2901

Figure 3. Sequence alignment of human PDE2A with PDE5A (pdb: 1XOZ.pdb). Stars refer to identical residues, whereas filled period refers to similar substitutions. R-Helices are indicated in underlined bold residues.

the side chains of Tyr827 and Phe830 residues. In addition, the molecular volume of Bay60-7550 structure is calculated to be ∼1240 Å3, larger than that of ND7001. The total volume of the PDE2AX-ray binding cavity is only about 1050 Å3. In comparison, the molecular volumes of the studied ligands were larger than that of the PDE2AX-ray active site cavity by ∼30 and 190 Å3 for ND7001 and Bay60-7550, respectively. Homology Modeling of Human PDE2A in a Bound State. The results discussed above reveal that the reported X-ray crystal structure21 of PDE2A in the unbound state is not suitable for molecular docking studies on the potent PDE2A inhibitors. Hence, we built a 3D model of PDE2A in a bound state (denoted by PDE2AHomology) through homology modeling using the highresolution X-ray crystal structure of PDE5A bound with a potent PDE5 inhibitor (Tadalafil)18 as the template. The amino acid sequence alignment revealed ∼32% identity between the human PDE2A and PDE5A (Figure 3). To investigate the quality of the refined PDE2AHomology structure, a variety of computational tests (described above) were performed. The first of these tests used Ramachandran plots to analyze the φ and ψ angle distributions for the refined model. The analysis produced by Procheck of the main-chain torsion angles of the PDE2AHomology revealed that a total of 91.0% of the residues are in the most favored regions of the Ramachandran plot with 8.3% in additionally allowed regions, giving a total of 99.3% in favorable or allowed regions. The results of the analysis suggest that the backbone conformation of our PDE2AHomology structure is reasonable. Other stereochemical parameters, such as the peptide planarity, bad nonbonded interactions, main-chain hydrogen bonding energy, and standard deviations of χ1 angle (i.e., the

first torsion angle of the side chain), were also examined. All these parameters obtained from the structural validation showed similarly good quality values for the PDE2AHomology structure compared to those obtained for the X-ray crystal structures of PDE2A and template PDE5 (Figure S1 in Supporting Information). The second test was to evaluate the quality of our PDE2AHomology structure through a detailed comparison of the packing environment of each residue in the modeled structure with the average packing environment of the same type of residue appeared in high-quality experimental structures deposited in the Protein Data Bank (PDB) using the Whatif program. A residue in a model structure with a score of -5.0 or worse usually indicates poor packing. The results of the Whatif quality report of the final PDE2AHomology, describing the evaluation of the structural integrity of the model, showed a z-score of -0.6. The z-scores calculated for individual residues of the PDE2AHomology and X-ray structures of PDE2A are provided as Supporting Information (Figure S2); the values all fall in the acceptable range for a valid structure (z-score >-5.0). Since there are even less residues with a score of -3.0 or lower in our PDE2AHomology model than that in the X-ray structure, the homology model has an acceptable packing quality in comparison with the available X-ray structure. The third test on the PDE2AHomology model was to apply energy criteria using the Prosa-II. We examined whether the interaction energy of each residue with the remainder of the protein is a negative value. Prosa-II energy plots for the PDE2AHomology model, X-ray structure of PDE2A, and the template X-ray structure of PDE5 are shown in Figure 4. It can

2902 J. Phys. Chem. B, Vol. 113, No. 9, 2009

Hamza and Zhan

Figure 4. (Top) Comparison of the Prosa-II quality report (z-score) of the PDE2AHomology 3D model with the PDE2A and PDE5A X-ray structures. (Bottom) Backbone atoms rms deviation of the whole PDE2AHomology protein and active site residues during the MD simulation.

be seen from Figure 4 that the overall results obtained for the X-ray structure of PDE2A and the PDE2AHomology structure are quite similar, as they have similar negative interaction energies. These results also suggest that the PDE2AHomology structure performs pretty much the same in the Prosa-II test as the solved X-ray structure of PDE2A. To validate the stability of the PDE2A structural model and to further test whether the PDE2AHomology structure modeled based on the high-resolution X-ray structure of human PDE5 is reasonable, we carried out a fully relaxed MD simulation on the PDE2AHomology structure in water. The MD simulation confirmed that the modeled PDE2AHomology structure was very stable (Figure 4), with a root-mean-square deviation (rmsd) of the backbone atoms remained almost intact during ∼3.5 ns of MD simulation. As expected, the PDE2AHomology model displayed a low rmsd fluctuation of ∼1.6 Å during the MD simulation. We also examined the rmsd of the backbone atoms of the residues in the PDE2AHomology active site and obtained a lower rmsd fluctuation during the MD simulation. The rmsd value converged to ∼1.2 Å after ∼1.7 ns, indicating that the active site structure of the protein model was very stable during the MD simulation. To sum up, the quality of the PDE2AHomology has been checked using four different types of criteria. The results show that the backbone conformation (Procheck), the residue interaction (Prosa-II), the residue contact (Whatif), and the stability of the protein folding (MD simulation) are all well within the ranges established for a reliable protein structure. Passing all tests on our PDE2AHomology structure suggests that we likely obtained an adequate 3D model of PDE2A to characterize its binding site and to explore PDE2A binding with the known inhibitors. Despite their similar overall quality, the difference between the refined PDE2AHomology structure and the X-ray structure is significant in the binding pocket. The structural superimposition between the 3D model and the X-ray structure revealed that the position of the LR14R15 loop and the R14 and R15 helices are more compact in the X-ray structure than that in the compu-

tational model. Indeed, a visual inspection of the PDE2AX-ray structure showed the presence of C-terminal residues (#902 to #915) in close contact with the solvent-accessible residue side chains of the LR14R15 loop and the R14 and R15 helices moiety. In addition, these residues, especially Leu907 and Leu913, are involved in hydrophobic interactions with Ile855 and Phe831 residues of the binding cavity, which induces a folding change of the LR14R15 loop such that it adopts an extended conformation. In contrast to the PDE2AHomology structure, the structural rearrangement of these residues in the PDE2AX-ray structure pulled inward Tyr827 and Phe830 side chains, giving rise to a more compact binding pocket. For another noticeable difference, as shown in Figure 5, the surface of the PDE2AHomology binding cavity is surrounded by two crevices (subsites) located at the interfaces between LR3R4 loop and R14 helix (upper domain) and between the R14 and R15 helices (lower domain), whereas in the PDE2AX-ray structure the subsite of the upper domain is fully filled by Tyr827 side chain. PDE2AHomology Binding with Bay60-7550 and ND7001. The computational experiment was designed to give an indication on the actual binding modes of PDE2A with Bay60-7550 and ND7001. The molecular docking and MD simulations were performed to reveal how these inhibitors bind with the protein and predict the corresponding PDE2AHomology-ligand binding free energies. To see how each of these inhibitors can anchor in the PDE2AHomology binding site pocket, we first inserted the inhibitor into the binding cavity and oriented in a favorable position by slightly rotating Phe830 and Y827 side chains and, thus, preventing any steric hindrance with the nearby residues of the pocket. After these minor structural modifications, the complex was submitted for the docking procedure as described above. Thus, the most favorable PDE2AHomology-Bay60-7550 binding mode was identified using clustering of the binding conformations (Figure 6). Molecular docking revealed few clusters of the binding modes, but only the primary cluster with the most populated conformation is associated with the expected strong binding between Phe830 and Ile826 side chains and this binding mode corresponds to the lowest interaction energy. The

Structure of PDE2 in a Bound State

J. Phys. Chem. B, Vol. 113, No. 9, 2009 2903

Figure 5. (Left) Most important residues of the PDE2A active site X-ray structure. The C-terminal residues are colored in brown. (Right) Superimposition of the PDE2AHomology model with the resolved X-ray structure (white color). Solvent-accessible surface (SAS) area of the active site is also displayed. The SAS area of the X-ray structure is represented in white meshed line.

Figure 6. (Left) View of the PDE2AHomology-Bay60-7550 complex with the ligand poses in the most populated cluster. The binding pocket of the receptor is represented by the SAS area. (Right) Average structure of the refined PDE2AHomology-Bay60-7550 complex. Bidentate hydrogen bonds between the ligand and the Gln859 side chain are displayed in yellow dashed line. (Bottom) Plots of MD-simulated internuclear distances, and rmsd for atomic positions of the ligand versus simulation time for PDE2AHomology binding with Bay60-7550. Traces D1 and D2 represent the bidentate H-bond distances between the oxygen and hydrogen of the ligand amide group and the γ-amide side chain of Gln859 residue. Traces D3 and D4 represent the distances between the C1 and C2 atoms (see Figure 1) of Bay60-7550 and the His656-Hε1 and Phe830-Cγ atoms. The graph represents the MD simulation with force constraint on the protein backbone.

results consistently show that the Bay60-7550 scaffold fits into the binding site cavity with the most conformations in close contact with Gln859 and Phe862 side chains of PDE2AHomology (Figure 6). So, Bay60-7550 binding with PDE2AHomology can be regarded as strong and specific as it reaches Gln859 side chain by hydrogen bonds. The PDE2AHomology-Bay60-7550 binding structure with the lowest interaction energy obtained from the docking was then refined by performing MD simulation. Inspection of the trajectory of the unconstrained MD simulation revealed that the positions of ligand heavy atoms remained stable after 1.1 ns of the MD simulation (see Figure S3 in Supporting Information). The jumps observed at ∼1 ns were due to the adoption of a

more stable conformation of Bay60-7550 in the pocket, producing a minor change on the flexible loops in the binding site. During the constrained MD simulation, the rmsd values of the ligand heavy atoms compared to the initial conformation were found to be stable and were rather small, i.e., 0.86 Å. The fluctuation of the two phenyl rings was described by the distances D3 and D4 in Figure 6. The average Bay60-7550 structure obtained from the MD simulation is shown in Figure 6, which reveals that Gln859 side chain and the carbonyl group of the ligand occupy the positions suitable for the formation of favorable hydrogen bonds. The amide moiety of the pyrazolopyrimidinone group of Bay607550 forms a bidentate hydrogen bond with the γ-amide group

2904 J. Phys. Chem. B, Vol. 113, No. 9, 2009

Hamza and Zhan

Figure 7. (Left) View of the PDE2AHomology-ND7001 complex with the ligand poses in the most populated cluster. The binding pocket of the receptor is represented by the SAS area. (Right) Average structure of the refined PDE2AHomology-ND7001 complex. Bidentate hydrogen bonds between the ligand and the Gln859 side chain are displayed in yellow dashed line. (Bottom) Plots of MD-simulated internuclear distances and rmsd for atomic positions of the ligand versus simulation time for PDE2AHomology binding with ND7001. Traces D1 and D2 represent the H-bond distances between the oxygen and hydrogen of the ligand amide group and the γ-amide side chain of Gln859 residue. Trace D3 represents the H-bond distance between the oxygen of the ligand methoxy group and the hydrogen amide side chain of Asn704 residue.

of Gln859, while the pyrazole ring forms hydrophobic interactions with the His656, Ile826, Leu858, Tyr827, and Phe830 side chains. We also note that, in the MD-simulated structure, the ligand was stabilized in the binding site by a network of π-π interactions (face to face) between Phe830 side chain and the two phenyl rings of the ligand and between and Phe862 side chain and the pyrazole scaffold. In addition, the nonsubstituted phenyl ring of the ligand was stabilized by a dipole-quadrupole interaction with His656 residue, while the metoxyphenyl ring is stabilized by a dipole-quadrupole interaction with Leu858 and Tyr827 residues. The similar molecular docking and MD simulation were also performed to examine PDE2AHomology binding with ND7001. The most frequent docking poses were found to cluster at hydrophobic regions and close to the conserved Gln859 residue (Figure 7). The obtained favorable binding structures with the best ligand fit were examined, and the structure with the lowest interaction energy was refined by MD simulation in water. Figure 7 shows that the MD trajectory of the simulated PDE2AHomology-ND7001 complex structure became stable after the first stage of MD simulation; the all-atom rms deviation of atomic positions of the ligand from the starting structure was only ∼0.5 Å. The stability of the PDE2AHomology-ND7001 binding structure was observed by checking some key internuclear distances between the ligand atoms and key amino acid residues of the protein. The side chain of the highly conserved Gln859 is the primary anchor point for the amide group of the ligand as reflected by the distances D1 and D2 in Figure 7, whose interactions are further stabilized by a hydrogen bond with the NH2 group of Asn704 side chain. These three hydrogen bonds were persistent during the MD simulation, with an average dynamic length close to ∼1.9 Å (which refers to the distances with the hydrogen in the hydrogen bonds).

Besides these important hydrogen bonds, several favorable interactions between the ligand and the enzyme are demonstrated in Figure 7. As observed in the PDE2AHomology-Bay60-7550 binding, the core structure of ND7001 is also lodged in the hydrophobic pocket involving Ile826 and Phe862, making a face-to-face π-π interaction with Phe862 side chain. Moreover, ND7001 is also stabilized by forming two tight dipole-quadrupole interactions between the methyl group of Ile826 side chain and the core of the ligand and between His656 side chain and the second phenyl group of the ligand. Further comparison of the simulated inhibitor-bound PDE2AHomology structures with the unbound form of the PDE2AXray led us to note two major structural differences between the inhibitor-bound and unbound states and obtain some insights into the possible structural motion of the enzyme active site. First, the binding of an inhibitor to the PDE2AHomology active site pocket is accompanied by a concerted closing motion of the pocket. A similar large motion can be seen in the comparison between the unbound form of PDE5 structure (2H40.pdb)63 and the bound form of several PDE5-ligand complexes (pdb codes: 2H44, 1XOZ, 1TBF, and 1RKB).63 The “induced-fit” movement in the PDE2A active site involves a loop consisting of residues 838 to 855 that acts as a lid over the inhibitor-binding site upon ligand binding. In addition, the core structure of the inhibitors shifted Phe830 side chain which was originally involved in covering the hydrophobic pocket of the lower binding domain (one of the subsites mentioned above) involving the R14 and R15 helices (Figure 5). The side chain of Phe830 rotates ∼90° from its original position toward the solvent in the unbound state. In the new position, Phe830 side chain is involved in π-π stacking interactions with the phenyl ring(s) of the ligand. The binding free energies for PDE2AHomology binding with Bay60-7550 and ND7001 determined by performing the MM-

Structure of PDE2 in a Bound State

J. Phys. Chem. B, Vol. 113, No. 9, 2009 2905

TABLE 1: Energetic Results (kcal/mol) Obtained from MM-PBSA Binding Free Energy Calculations on PDE2A Binding with Bay60-7550 and ND7001 PDE2AHomology PDE2AX-rayfbound state PDE2AHomology -Bay60-7550 - Bay60-7550 -ND7001 ele ∆Eint vdw ∆Eint ∆EMM ∆Esol ∆Ebind -T∆S ∆Gbind(calc) ∆Gbind(expt)a

-18.2 (1.4) -66.0 (1.2) -84.2 (1.3) 44.8 (1.7) -39.4 (1.4) 26.5 (1.1) -12.9 (1.2)

-19.9 (1.2) -66.4 (1.4) -86.3 (1.3) 47.1 (1.2) -39.2 (1.6) 26.5 (1.2) -12.7 (1.4) -11.7

-29.2 (0.9) -45.5 (0.7) -74.7 (0.7) 46.9 (0.8) -27.8 (0.8) -18.5 (1.4) -9.3 (1.1) -10.4

a

Experimental values derived from the available IC50 values in refs 23 and 25.

PBSA calculations are summarized in Table 1. The calculated binding free energies are -12.9 kcal/mol for Bay60-7550 and -9.3 kcal/mol for ND7001. The calculated binding free energies are close to the corresponding experimental data, i.e., -11.7 kcal/mol for Bay60-7550 and -10.4 kcal/mol for ND7001, suggesting that the MM-PBSA binding free energy calculations are reasonable for these protein-ligand binding systems. In consideration of the binding free energy components, we note that the overall van der Waals interaction with Bay60-7550 (-66.0 kcal/mol) is stronger than that with ND7001 (-45.5 kcal/mol) by ∼21 kcal/mol, whereas the overall electrostatic interaction with Bay60-7550 (-18.2 kcal/mol) is weaker than that with ND7001 (-29.2 kcal/mol) by ∼11 kcal/mol. Indeed, the PDE2A-ND7001 binding mode is due to strong electrostatic interactions (hydrogen bonds) with Asn704 and Gln859 side chains. The strong van der Waals and electrostatic interactions are partly compensated by the unfavorable solvation energy term. The calculated energetic results suggest that it is crucial for rational design of a high-affinity inhibitor of PDE2A to achieve an optimal van der Waals interaction between the ligand and the protein active site and to suffer less desolvation penalty. Further Examination of the Inhibitor-Bound Form of PDE2A Structure. Taking into account all of the abovediscussed structures concerning the bound and unbound states of PDE2A and their binding with inhibitors, we wanted to know whether we could also model the 3D structure of the inhibitorbound state of PDE2A starting from the X-ray crystal structure, i.e., PDE2AX-ray, associated with the unbound state in which the binding pocket was in the collapsed form. As discussed above, in the PDE2AHomology structure, the active site is extended to include the aforementioned subsite under loop LR14R15 (residues 838 to 855) compared to the PDE2AX-ray structure. For convenience, the 3D model of the bound state of PDE2A modeled from the PDE2AX-ray structure is denoted by PDE2AXX-rayfbound state structure, we first rayfbound state. To model the PDE2A rebuilt LR14R15 loop of the PDE2AX-ray by using the structure of LR14R15 loop in the X-ray crystal structure of PDE5 as a template. The obtained initial structure was then relaxed by carrying out MD simulation in water. Taking into account the protein flexibility, we may consider PDE2AX-rayfbound state to be an ensemble of many protein structures.64,65 The structural ensemble was obtained from the snapshots of the last 1.2 ns of the MD trajectory. Since it was observed from the MD simulation that the CR rmsd of the LR14R15 was found to be less than ∼0.9 Å, we decided to average those conformations into the final conformation of the loop. Figure 8 shows the superimposition of six representative snapshot structures of PDE2AX-rayfbound state. Thus during the

docking, the ligand interacted with a mean-field due to the average of protein structures. The ensemble structures obtained from the MD simulation demonstrated a stabilized LR14R15 loop in an extended conformation. During the MD simulation, Phe830, Tyr827, and Phe831 side chains rotated automatically toward outside of the active site cavity. As mentioned above, in the PDE2AHomology structure, the solvent-accessible surface area of PDE2AX-rayfbound state displays a larger binding cavity compared to the X-ray structure, and the active site cavity in the PDE2AHomology structure includes additional two deep subsites that may be filled by the two phenyl moieties of Bay60-7550 (Figure 8). To construct a reasonable initial PDE2AX-rayfbound state-Bay60-7550 complex structure, the refined PDE2AHomology-Bay60-7550 complex structure discussed above was first superimposed to the PDE2AX-rayfbound state structure and then the Bay60-7550 molecule from the complex was assembled to the PDE2AX-rayfbound state structure. The obtained initial structure of PDE2AX-rayfbound state-Bay60-7550 complex was reaccommodated by slightly rotating Ile826, Tyr827, and Phe830 side chains and was energy-minimized. The complex structure was then refined during ∼2.5 ns of a fully relaxed MD simulation in water and the MD simulation continued for 1.3 ns with force constrain (50 kcal mol-1 Å-2) on backbone atoms. The results obtained from the MD simulation on the complex in water reveal that the protein structure in the MD-simulated PDE2AX-rayfbound state-Bay60-7550 complex is very close to that in the MD-simulated PDE2AHomology-Bay60-7550 complex. Especially, the simulated binding mode of Bay60-7550 in the PDE2AX-rayfbound state structure was very similar to that in the PDE2AHomology model (see Figure S4 in Supporting Information). For example, the strongest interactions presented by the pyrazolopyrimidinone group with Gln859 and Phe862 residues in the PDE2AHomology-Bay60-7550 complex structure were well reflected in the PDE2AX-rayfbound state-Bay60-7550 complex structure. Hence, it is not surprising to find out that the MMPBSA-calculated binding free energy (-12.7 kcal/mol) for Bay60-7550 binding with PDE2AX-rayfbound state is also close to that (-12.9 kcal/mol) with PDE2AHomology (see Table 1). Both computational values (-12.7 and -12.9 kcal/mol) are reasonably close to the experimentally derived value of -11.7 kcal/ mol. All of the computational results suggest that the bound form of the PDE2A structure not only can be modeled from the corresponding bound form of PDE5 structure but also can be transformed from the unbound form of the PDE2A structure by adjusting the LR14R15 loop conformation and relaxing the entire structure. The structural change from the unbound to the bound state leads to a substantial variation in the size of the pocket but does not change the general structural feature of the catalytic site, as shown in recently reported crystal structures of PDE5 complexed with ligands in different sizes.63 The similar phenomena was observed in the structural transition from the apo to holo forms of proteins and/or the ligand-dependent structural changes involved in most protein kinases,66-68 in transcriptional activation of the nuclear hormone receptor (NR) family69,70 and in estrogen receptor-R.71,72 The ability of PDE2A to adopt different structural forms highlights the extreme flexibility and adaptability of the catalytic domain of the enzyme, allowing the molding of the protein around its ligand and consequently the formation of the proper ligand-binding pocket. As suggested by the significantly higher B-factors in the LR14R15 “lid” loop of the resolved X-ray crystal structure, the unliganded PDE2A protein is expected to open

2906 J. Phys. Chem. B, Vol. 113, No. 9, 2009

Hamza and Zhan

Figure 8. (Left) Superimposition of the structural ensemble of the PDE2AX-rayfbound state snapshots (green) obtained during the MD simulation with the resolved PDE2AX-ray structure (purple). (Right) Superimposition of the average of the PDE2AX-rayfbound state structural snapshots (green) with the PDE2AX-ray structure (purple). Solvent-accessible surface (SAS) area of the active site is also displayed. (Bottom) Plots of MD-simulated internuclear distances, and rmsd for atomic positions of the ligand vs simulation time for PDE2AX-rayfbound state binding with Bay60-7550. Traces D1 and D2 represent the H-bond distances between the oxygen and hydrogen of the ligand amide group and the γ-amide side chain of Gln859 residue. Traces D3 and D4 represent the distances between the C1 and C2 atoms (see Figure 1) of Bay60-7550 and the His656-H1 and Phe830-Cγ atoms. The graph represents the MD simulation with force constraint on the protein backbone.

and close frequently in solution. When a substrate or inhibitor binds to the binding domain of the open state, it will be trapped when the protein closes. In addition, all of the results suggest an interesting new concept, i.e., the ligand-dependent binding pocket, in the field of PDE proteins. In fact, although we propose the concept of ligand-dependent binding pocket of a PDE protein for the first time in this study, it actually can be seen from the recently resolved X-ray crystal structures of other PDE proteins18,20,63 that certain local adaptations of the catalytic active site of PDE to its ligand can take place. These structural adaptations/modifications of PDE proteins might lead to a substantial variation in the size of the binding pocket, but do not the general feature of the catalytic site structure, as shown in the crystal structures of PDE4/5 complexed with ligands with different sizes.18,20,63

obtained in this study are very close to each other and have led to very similar PDE2A-inhibitor binding modes and binding free energies. All of the calculated binding free energies are in good agreement with available experimental data. The ability of PDE2A to adopt different structural forms while maintaining the general structural feature of the catalytic site suggests that the binding pocket of PDE2A is very flexible and may accommodate a variety of ligands with distinctly different molecular structures, which leads us to better understand why PDE2 can catalyze the hydrolysis of both cAMP and cGMP whereas PDE4 and PDE5 can only respectively catalyze the hydrolysis of cAMP and cGMP. The general structural insights, PDE2A model of PDE2A in the bound state, and detailed structures of the PDE2A-inhibitor binding structures obtained in this study will be valuable for future rational design of novel, potent inhibitors of PDE2A as therapeutic agents.

Conclusion Extensive molecular modeling, docking, and dynamics simulations have demonstrated that the structure of PDE2A catalytic domain may exist in two different forms corresponding to the inhibitor-bound and unbound states of the enzyme. The most remarkable structural difference between the two states is associated with the LR14R15 loop which exists in a more extended conformation in the bound state compared to the unbound state. The structural change from the unbound state to the bound state leads to a substantial variation in the size of the binding pocket, but does not affect the general structural feature of the catalytic site. The available X-ray crystal structure of PDE2A was in the unbound state and, thus, is not suitable for molecular docking studies on PDE2A binding with its known potent inhibitors. The PDE2A structure of the bound state can be modeled starting from either the X-ray crystal structure of PDE5 with an inhibitor bound in the binding pocket or the available X-ray crystal structure of PDE2A in the unbound state. The two 3D models

Acknowledgment. This research was supported in part by Kentucky Science and Engineering Foundation (grant KSEF925-RDE-008). The authors also acknowledge the Center for Computational Sciences (CCS) at University of Kentucky for supercomputing time on IBM X-series Cluster with 1360 processors. Supporting Information Available: Six figures concerning more detailed information about the homology modeling, MDsimulated structures, and MD trajectories of PDE2A binding with the inhibitors. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Burns, F.; Zhao, A.; Beavo, J. A. Cyclic nucleotide phosphodiesterases: gene complexity, regulation by phosphorylation and physiological implications. AdV. Pharmacol. 1996, 36, 39–48.

Structure of PDE2 in a Bound State (2) Beavo, J.; Conti, M.; Heaslip, R. J. Multiple cyclic nucleotide phosphodiesterase. Mol. Pharmacol. 1994, 46, 399–405. (3) Conti, M.; Jin, S. L. The molecular biology of cyclic nucleotide phosphodiesterases. Prog. Nucleic Acid Res. Mol. Biol. 1999, 63, 1–38. (4) Houslay, M. D.; Schafer, P.; Zhang, K. Y. Keynote review: phosphodiesterase-4 as a therapeutic target. Drug DiscoV. Today 2005, 10, 1503–1519. (5) Georget, M.; Mateo, P.; Vandecasteele, G.; Lipskaia, L.; Defer, N.; Hanoune, J.; Hoerter, J.; Lugnier, C.; Fischmeister, P. Cyclic AMP compartmentation due to increased cAMP-phosphodiesterase expression in transgenic mice with a cardiac-directed expression of the human adenylyl cyclase type 8 (AC8). FASEB J. 2003, 17, 1380–1391. (6) Karpen, J. W.; Rich, T. C. High-resolution measurements of cyclic adenosine monophosphate signals in 3D domains. Methods Mol. Biol. 2005, 307, 15–26. (7) Loughney, K.; Ferguson, K. Identification and quantification of PDE isozymes and subtypes by molecular biological methods. In Phosphodiesterase Inhibitors; Schudt, C., Dent, G., Rabe, K. F., Eds.; Academic Press: San Diego, CA, 1996; pp 1-19. (8) Conti, M.; Jin, S. L. The molecular biology of cyclic nucleotide phosphodiesterases. Prog. Nucleic Acid Res. Mol. Biol. 1999, 63, 1–38. (9) Soderling, S. H.; Beavo, J. A. Regulation of cAMP and cGMP signaling: new phosphodiesterases and new functions. Curr. Opin. Cell. Biol. 2000, 12, 174–179. (10) Rosman, G. J.; Martins, T. J.; Sonnenburg, W. K.; Beavo, J. A.; Ferguson, K.; Loughney, K. Isolation and characterization of human cDNAs encoding a cGMP-stimulated 3′,5′-cyclic nucleotide phosphodiesterase. Gene 1997, 191, 89–95. (11) Rivet-Bastide, M.; Vandecasteele, G.; Hatem, S.; Verde, I.; Be´nardeau, A.; Mercadier, J. J.; Fischmeister, R. cGMP-stimulated cyclic nucleotide phosphodiesterase regulates the basal calcium current in human atrial myocytes. J. Clin. InVest. 1997, 99, 2710–2718. (12) Sadhu, K.; Hensley, K.; Florio, V. A.; Wolda, S. L. Differential expression of the cyclic GMP-stimulated phosphodiesterase PDE2A in human venous and capillary endothelial cells. J. Histochem. Cytochem. 1999, 47, 895–906. (13) Martinez, S. E.; Wu, A. Y.; Glavas, N. A.; Tang, X. B.; Turley, S.; Hol, W. G. J.; Beavo, J. A. The two GAF domains in phosphodiesterase 2A have distinct roles in dimerization and in cGMP binding. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 13260–13265. (14) Xu, R. X.; Rocque, W. J.; Lambert, M. H.; Vanderwall, D. E.; Luther, M. A.; Nolte, R. T. Crystal structure of the catalytic domain of phospohodiesterase 4B complexed with AMP, 8-Br-AMP and rolipram. J. Mol. Biol. 2004, 337, 355–365. (15) Houslay, M. D.; Schafer, P.; Zhang, K. Y. Keynote review: phosphodiesterase-4 as a therapeutic target. Drug DiscoV. Today 2005, 10, 1503–1519. (16) Raja, S. G.; Nayak, S. H. Sildenafil: emerging cardiovascular indications. Ann. Thorac. Surg. 2004, 78, 1496–1506. (17) Schror, K. The pharmacology of cilostazol. Diabetes Obes. Metab. 2002, 4 (Suppl 2), S14-S19. (18) Card, G. L.; England, B. P.; Suzuki, Y.; Fong, D.; Powell, B.; Lee, B.; Luu, C.; Tabrizizad, M.; Gillette, S.; Ibrahim, P. N.; Artis, D. R.; Gideon Bollag, G.; Milburn, M. V.; Kim, S. H.; Schlessinger, J.; Zhang, K. Y. J. Structural basis for the activity of drugs that inhibit phosphodiesterases. Structure 2004, 12, 2233–2247. (19) Zhang, K. Y. J.; Card, G. L.; Suzuki, Y.; Artis, D. R.; Fong, D.; Gillette, S.; Hsieh, D.; Neiman, J.; West, B. L.; Zhang, C.; Milburn, M. V.; Kim, S.-H.; Schlessinger, J.; Bollag, G. A glutamine switch mechanism for nucleotide selectivity by phosphodiesterases. Mol. Cell 2004, 15, 279– 286. (20) Jeon, Y. H.; Heo, Y. S.; Kim, C. M.; Hyun, Y. L.; Lee, T. G.; Ro, S.; Cho, J. M. Phosphodiesterase: overview of protein structures, potential therapeutic applications and recent progress in drug development. Cell. Mol. Life Sci. 2005, 62, 1198–1220. (21) Iffland, A.; Kohls, D.; Low, S.; Luan, J.; Zhang, Y.; Kothe, M.; Cao, Q.; Kamath, A. V.; Ding, Y. H.; Ellenberger, T. Structural determinants for inhibitor specificity and selectivity in PDE2A using the wheat germ in vitro translation system. Biochemistry 2005, 23, 8312. (22) Beavo, J. A. Cyclic nucleotide phosphodiesterases: functional implications of multiple isoforms. Physiol. ReV. 1995, 75, 725–48. (23) Boess, F. G.; Hendrix, M.; van der Staay, F. J.; Erb, C.; Schreiber, R.; van Staveren, W.; de Vente, J.; Prickaerts, J.; Blokland, A.; Koenig, G. Inhibition of phosphodiesterase 2 increases cGMP neuronal plasticity and memory performance. Neuropharmacology 2004, 47, 1081–92. (24) Podzuweit, T.; Nennstiel, P.; Muller, A. Isozyme selective inhibition of cGMP-stimulated cyclic nucleotide phosphodiesterases by erythro-9-(2hydroxy-3-nonyl) adenine. Cell Signal 1995, 7, 733–8. (25) Abarghaz, M.; Biondi, S.; Duranton, J.; Limanton, E.; Mondadori, C. Benzodiazepine derivatives, their preparation and the therapeutic use thereof. European Patent EP1749824, 2005. (26) Benson, D. A.; Karsch-Mizrachi, I.; Lipman, D. J.; Ostell, J.; Wheeler, D. L. Genbank. Nucleic Acids Res. 2000, 28, 15–8.

J. Phys. Chem. B, Vol. 113, No. 9, 2009 2907 (27) Bernstein, F. C.; Koetzle, T. F.; Williams, G. J.; Meyer, E. E., Jr.; Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. The Protein Data Bank: a computer-based archival file for macromolecular structures. J. Mol. Biol. 1977, 112, 535–42. (28) Altschul, S. F.; Gish, W.; Miller, W.; Myers, E. W.; Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 1990, 215, 403–10. (29) Fiser, A.; Do, R. K.; Sali, A. Modeling of loops in protein structures. Protein Sci. 2000, 9, 1753–1773. (30) Marti-Renom, M. A.; Stuart, A. C.; Fiser, A.; Sanchez, R.; Melo, F.; Sali, A. Comparative protein structure modeling of genes and genomes. Annu. ReV. Biophys. Biomol. Struct. 2000, 29, 291–325. (31) Sali, A.; Blundell, T. L. Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 1993, 234, 779–815. (32) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A program for macromolecular energy minimization and dynamics calculations. J. Comput. Chem. 1983, 4, 187–217. (33) Sutton, R. B.; Rasshauer, D.; Jahn, R.; Brunger, A. T. Crystal structure of a SNARE complex involved in synaptic exocytosis at 2.4 A resolution. Nature 1998, 395, 347–353. (34) Flohil, J. A.; Vriend, G; Berendsen, H. J. C. Completion and refinement of 3-D homology models with restricted molecular dynamics: Application to targets 47, 58, and 111 in the CASP modeling competition and posterior analysis. Proteins-Struct., Funct., Genet. 2002, 48, 593– 604. (35) Laskowski, R. A.; MacArthur, M. W.; Moss, D. S.; Thornton, J. M. PROCHECK: a program to check the stereochemical quality of protein structures. J. Appl. Crystallogr. 1993, 26, 283–291. (36) Morris, A. L.; MacArthur, M. W.; Hutchinson, E. G.; Thornton, J. M. Stereochemical quality of protein structure coordinates. Proteins 1992, 12, 345–364. (37) Sippl, M. J. Recognition of errors in three-dimensional structures of proteins. Proteins 1993, 17, 355–362. (38) Vriend, G. WHAT IF: a molecular modeling and drug design program. J. Mol. Graph. 1990, 8, 52. (39) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millan, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malich, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzales, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andreas, J. L.; HeadGordon, M.; Reploge, E. S.; Pople, J. A. Gaussian 03, reVision B-03; Gaussian, Inc.: Pittsburgh, PA, 2003. (40) (a) Cieplak, P.; Cornell, W. D.; Bayly, C. I.; Kollman, P. A. Application of the multimolecule and multiconformational RESP methodology to biopolymers: charge derivation for DNA, RNA, and proteins. J. Comput. Chem. 1995, 16, 1357. (b) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys. Chem. 1993, 97, 10269. (41) Case, D. A.; Darden, T. A.; Cheatham III, T. E.; Simmerling, C. L.; Wang, J.; Duke, et al. AMBER 9; University of California: San Francisco, 2006. (42) Hamza, A.; Zhan, C-G. How can (-)-epigallocatechin gallate from green tea prevent HIV-1 infection? mechanistic insights from computational modeling and the implication for rational design of anti-HIV-1 entry inhibitors. J. Phys. Chem. B 2006, 110, 2910–2917. (43) Bargagna-Mohan, P.; Hamza, A.; Kim, Y. E.; Khuan, Y.; Ho, A.; Mor-Vaknin, N.; Wendschlag, N.; Liu, J.; Evans, R. M.; Markovitz, D. M.; Zhan, C-G.; Kim, K. B.; Mohan, R. The tumor inhibitor and antiangiogenic agent withaferin A targets the intermediate filament protein vimentin. Chem. Biol. 2007, 14]?>, 623–34. (44) Zhang, T.; Hamza, A.; Cao, X.; Wang, B.; Yu, S.; Zhan, C. G.; Sun, D. A novel Hsp90 inhibitor to disrupt Hsp90/Cdc37 complex against pancreatic cancer cells. Mol. Cancer Ther. 2008, 7, 162–70. (45) Kuntz, I. D.; Meng, E. C.; Shoichet, B. K. Structure-based molecular design. Acc. Chem. Res. 1994, 27, 117. (46) Pearlman, D. A., et al. AMBER 4.1; University of California: San Francisco, CA, 1995. (47) Harvey, S. C. Treatment of electrostatic effects in macromolecular modeling. Proteins 1989, 5, 78–92. (48) Guenot, J.; Kollman, P. A. Molecular dynamics studies of a DNAbinding protein: 2. An evaluation of implicit and explicit solvent models for the molecular dynamics simulation of the Escherichia coli trp repressor. Protein Sci. 1992, 1, 1185–1205. (49) Gatchell, D. W.; Dennis, S.; Vajda, S. Discrimination of near-native protein structures from misfolded models by empirical free energy functions. Proteins 2000, 41, 518–534.

2908 J. Phys. Chem. B, Vol. 113, No. 9, 2009 (50) Kozakov, D.; Clodfelter, K. H.; Vajda, S.; Camacho, C. J. Optimal clustering for detecting near-native conformations in protein docking. Biophys. J. 2005, 89, 867–875. (51) Jorgensen, W. L. Quantum and statistical mechanical studies of liquids. 10. Transferable intermolecular potential functions for water, alcohols, and ethers. Application to liquid water. J. Am. Chem. Soc. 1981, 103, 335. (52) Pan, Y.; Gao, D.; Yang, W.; Cho, H.; Yang, G-F.; Tai, H-H.; Zhan, C-G. Computational redesign of human butyrylcholinesterase for anticocaine medication. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 16656. (53) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: an N · log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089. (54) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684. (55) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327. (56) Abdulhameed, M. D. M.; Hamza, A.; Zhan, C-G. Microscopic modes and free energies of 3-phosphoinositide dependent kinase-1 (PDK1) binding with celecoxib and other inhibitors. J. Phys. Chem. B. 2006, 110, 26365–26374. (57) Kollman, P.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.; Srnivasan, J.; Case, D.; Cheatham, T. E. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889–897. (58) Hou, T.; Guo, S.; Xu, X. Predictions of binding of a diverse set of ligands to Gelatinase-A by a combination of molecular dynamics and continuum solvent models. J. Phys. Chem. B 2002, 106, 5527–5535. (59) (a) Xiong, Y.; Li, Y.; He, H; Zhan, C. G. Theoretical calculation of the binding free energies for pyruvate dehydrogenase E1 binding with ligands. Bioorg. Med. Chem. Lett. 2007, 17, 5186–5190. (b) Pan, Y.; Gao, D.; Zhan, C.-G. Modeling the catalysis of anti-cocaine catalytic antibody: Competing reaction pathways and free energy barriers. J. Am. Chem. Soc. 2008, 130, 5140–5149. (60) Xu, R. X.; Hassell, A. M.; Vanderwall, D.; Lambert, M. H.; Holmes, W. D.; Luther, M. A.; Rocque, W. J.; Milburn, M. V.; Zhao, Y.; Ke, H.; Nolte, R. T. Atomic structure of PDE4: Insights into phosphodiesterase mechanism and specificity. Science 2000, 288, 1822–1825.

Hamza and Zhan (61) Lee, M. E.; Markowitz, J.; Lee, J.-O.; Lee, H. Crystal structure of phosphodiesterase 4D and inhibitor complex. FEBS Lett. 2002, 530, 58– 53. (62) Sung, B. J.; Hwang, K. Y.; Jeon, Y. H.; Lee, J. I.; Heo, Y. S.; Kim, J. H.; Moon, J.; Yoon, J. M.; Hyun, Y. L.; Kim, E.; Eum, S. J.; Park, S. Y.; Lee, J. O.; Lee, T. G.; Ro, S.; Cho, J. M. Structure of the catalytic domain of human phosphodiesterase 5 with bound drug molecules. Nature 2003, 425, 98–102. (63) Wang, H.; Liu, Y.; Huai, Q.; Cai, J.; Zoraghi, R.; Francis, S. H.; Corbin, J. D.; Robinson, H.; Xin, Z.; Lin, G.; Ke, H. Multiple conformations of phosphodiesterase-5: implications for enzyme function and drug development. J. Biol. Chem. 2006, 281, 21469–79. (64) Knegtel, R. M. A.; Kuntz, I. D.; Oshiro, C. M. Molecular docking to ensembles of protein structures. J. Mol. Biol. 1997, 266, 424–440. (65) Osterberg, G. M.; Morris, M. F.; Sanner, A. J.; Goodsell, O.; Goodsell, D. S. Automated docking to multiple target structures: incorporation of protein mobility and structural water heterogeneity in AutoDock. Proteins 2002, 46, 34–40. (66) Frost, J. A.; Khokhlatchev, A.; Stippec, S.; White, M. A.; Cobb, M. H. Differential effects of PAK1-activating mutations reveal activitydependent and -independent effects on cytoskeletal regulation. J. Biol. Chem. 1998, 273, 28191–28198. (67) Lei, M.; Lu, W.; Meng, W.; Parrini, M. C.; Eck, M. J.; Mayer, B. J.; Harrison, S. C. Structure of PAK1 in an autoinhibited conformation reveals a multistage activation switch. Cell 2000, 102, 387–397. (68) Hoffman, G. R.; Cerione, R. A. Flipping the switch: the structural basis for signaling through the CRIB motif. Cell 2000, 102, 403–406. (69) Moras, D.; Gronemeyer, H. Review: The nuclear receptor ligandbinding domain: structure and function. Curr. Opin. Cell. Biol. 1998, 10, 384–91. (70) Klaholz, B. P.; Renaud, J. P.; Mitschler, A.; Zusi, C.; Chambon, P.; Gronemeyer, H.; Moras, D. Conformational adaptation of agonists to the human nuclear receptor RAR gamma. Nat. Struct. Biol. 1998, 5, 199– 202. (71) Brzozowski, A. M.; Pike, A. C.; Dauter, Z.; Hubbard, R. E.; Bonn, T.; Engstro¨m, O.; Ohman, L.; Greene, G. L.; Gustafsson, J. A.; Carlquist, M. Molecular basis of agonism and antagonism in the oestrogen receptor. Nature 1997, 389, 753–8. (72) Shiau, A. K.; Barstad, D.; Loria, P. M.; Cheng, L.; Kushner, P. J.; Agard, D. A.; Greene, G. L. The structural basis of estrogen receptor/ coactivator recognition and the antagonism of this interaction by tamoxifen. Cell 1998, 95, 927–37.

JP8082612