Biophysical Joumal Volume 93 September 2007 1431-1441 143 Dopamine D1 Receptor Agonist and D2 Receptor Antagonist Effects of the Natural Product (-Stepholidine: Molecular Modeling and dynamics simulations Wei Fu, t Jianhua Shen, *Xiaomin Luo, *Weiliang Zhu, Jiagao Cheng, *Kunqian Yu, James M. Briggs, Guozhang Jin, * Kaixian Chen, and Hualiang Jiang Discovery and Design Centre, State Key stitute of Materia Medica, Shangh titutes for cal Sciences, Chinese Academy of Scie nghai 201203, People' blic of China; Department of Medicinal Chemistry of Pharmacy, Fudan University, Shang Department of Biology and Biochemistry, Universit on, Houston, Texas, 77204-5001; and sSchool of Pharmacy, East China University of Science and Technology, Shanghai 200237, People's Republic of China ABSTRACT (-)Stepholidine(SPD), an active ingredient of the Chinese herb Stephania, is the first compound found to have dual function as a dopamine receptor D1 agonist and D2 antagonist. Insights into dynamical behaviors of D1 and D2 receptors and their interaction modes with SPD are crucial in understanding the structural and functional characteristics of dopamine receptors. In this study a computational approach, integrating protein structure prediction, automated molecular docking and molecular dynamics simulations were employed to investigate the dual action mechanism of SPd on the D1 and D2 receptors, with the eventual aim to develop new drugs for treating diseases affecting the central nervous system such as schizophrenia The dynamics simulations revealed the surface features of the electrostatic potentials and the conformational"open-closed process of the binding entrances of two dopamine receptors. Potential binding conformations of D1 and D2 receptors were obtained, and the D1-SPD and D2-SPD complexes were generated, which are in good agreement with most of experimental data. The D1-SPD structure shows that the K-167_EL-2-E-302_EL-3(EL-2: extracellular loop 2; EL-3: extracellular loop 3 )salt bridge plays an important role for both the conformational change of the extracellular domain and the binding of SPD. Based on our modeling and simulations, we proposed a mechanism of the dual action of SPD and a subsequent signal transduction model. Further mutagenesis and biophysical experiments are needed to test and improve our proposed dual action mechanism of SPD and signal transduction model. INTRODUCTION During the last decade, numerous efforts have been under- mental limitations and the lack of three-dimensional (3D) taken to discover drugs for treating psychomotor diseases, structures including the debilitating mental illness schizophrenia, Antagonists of the D2 receptor are believed to be potential which affects 0.5-1.5% of the worldwide population(1-4). drugs against the psychomotor diseases(9). Unfortunately, has been established that dopamine receptors(DRs)are these antagonists (e. g, the neuroleptic haloperidol) have primary targets for developing drugs to treat these diseases. severe mechanism-related side effects, including induction DRs belong to the G-protein coupled receptor(GPCR) of acute extrapyramidal symptoms(EPS), tardive dyskinesia, superfamily, transferring signals into cells through guanine and problems of galactorrhea due to increase in prolactin nucleotide-binding regulatory G-proteins(4). The DRs can release (10). In contrast, atypical antipsychotics with less be classified into two major subfamilies, DI and D2 re- EPS are effective in patients who are unresponsive to ceptors, according to their structural, pharmacological, and classical agents and may also have advantages in treating the functional characteristics (4). Their functional domains more resistant negative symptoms of schizophrenia(10). The have been defined, and the binding-site crevices have been use of atypical antipsychotics was, however, found to have a identified by the substituted-cysteine-accessibility method high occurrence of a potentially fatal blood disorder called (5-7). Although some studies revealed the structural fea- agranulocytosis(10). Therefore, discovering effective, safe tures of DI and D2 receptors that underlie their partic- antipsychotics remains a high priority(4, 9, 10) ular biophysical and pharmacological properties (6, 8). The pathogenesis of schizophrenia was suggested to be many problems still remain unresolved due to experi- related to dysfunction of the DI receptor in the medial pre- frontal cortex(mPFC), which is accompanied by D2 receptor hyperactivity in subcortical regions such as the ventral te ch20.2007 mental area (VTA)and the nucleus accumbens(NAc)(4) Address reprint requests to Prof. Hualiang Weiliang Zhu, The DI dysfunction was suggested to be responsible for the Shanghai Institute of Materia Medica. Chinese Jiang and of Sciences, 555 negative symptoms of schizophrenia, whereas the hyperacti hong zhi road, Zhangjiang Hi-Tech Park, Shanghai 201203 People' s ity of the D2 receptor might lead to the positive symptoms of hejiang@mail shcc.ac cn this disorder (4,5, 11). Based on this hypothesis, antipsychotic o 2007 by the Biophysical Society 0006-34950709/1431/11S200 doi:10.1529 biophys106.088500
Dopamine D1 Receptor Agonist and D2 Receptor Antagonist Effects of the Natural Product (2)–Stepholidine: Molecular Modeling and Dynamics Simulations Wei Fu,*y Jianhua Shen,* Xiaomin Luo,* Weiliang Zhu,* Jiagao Cheng,* Kunqian Yu,* James M. Briggs,z Guozhang Jin,* Kaixian Chen,* and Hualiang Jiang*§ *Drug Discovery and Design Centre, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Shanghai Institutes for Biological Sciences, Chinese Academy of Sciences, Shanghai 201203, People’s Republic of China; y Department of Medicinal Chemistry, School of Pharmacy, Fudan University, Shanghai 200032, People’ Republic of China; z Department of Biology and Biochemistry, University of Houston, Houston, Texas, 77204-5001; and § School of Pharmacy, East China University of Science and Technology, Shanghai 200237, People’s Republic of China ABSTRACT ()–Stepholidine (SPD), an active ingredient of the Chinese herb Stephania, is the first compound found to have dual function as a dopamine receptor D1 agonist and D2 antagonist. Insights into dynamical behaviors of D1 and D2 receptors and their interaction modes with SPD are crucial in understanding the structural and functional characteristics of dopamine receptors. In this study a computational approach, integrating protein structure prediction, automated molecular docking, and molecular dynamics simulations were employed to investigate the dual action mechanism of SPD on the D1 and D2 receptors, with the eventual aim to develop new drugs for treating diseases affecting the central nervous system such as schizophrenia. The dynamics simulations revealed the surface features of the electrostatic potentials and the conformational ‘‘open-closed’’ process of the binding entrances of two dopamine receptors. Potential binding conformations of D1 and D2 receptors were obtained, and the D1-SPD and D2-SPD complexes were generated, which are in good agreement with most of experimental data. The D1-SPD structure shows that the K-167_EL-2-E-302_EL-3 (EL-2: extracellular loop 2; EL-3: extracellular loop 3) salt bridge plays an important role for both the conformational change of the extracellular domain and the binding of SPD. Based on our modeling and simulations, we proposed a mechanism of the dual action of SPD and a subsequent signal transduction model. Further mutagenesis and biophysical experiments are needed to test and improve our proposed dual action mechanism of SPD and signal transduction model. INTRODUCTION During the last decade, numerous efforts have been undertaken to discover drugs for treating psychomotor diseases, including the debilitating mental illness schizophrenia, which affects 0.5–1.5% of the worldwide population (1–4). It has been established that dopamine receptors (DRs) are primary targets for developing drugs to treat these diseases. DRs belong to the G-protein coupled receptor (GPCR) superfamily, transferring signals into cells through guanine nucleotide-binding regulatory G-proteins (4). The DRs can be classified into two major subfamilies, D1 and D2 receptors, according to their structural, pharmacological, and functional characteristics (4). Their functional domains have been defined, and the binding-site crevices have been identified by the substituted-cysteine-accessibility method (5–7). Although some studies revealed the structural features of D1 and D2 receptors that underlie their particular biophysical and pharmacological properties (6,8), many problems still remain unresolved due to experimental limitations and the lack of three-dimensional (3D) structures. Antagonists of the D2 receptor are believed to be potential drugs against the psychomotor diseases (9). Unfortunately, these antagonists (e.g., the neuroleptic haloperidol) have severe mechanism-related side effects, including induction of acute extrapyramidal symptoms (EPS), tardive dyskinesia, and problems of galactorrhea due to increase in prolactin release (10). In contrast, atypical antipsychotics with less EPS are effective in patients who are unresponsive to classical agents and may also have advantages in treating the more resistant negative symptoms of schizophrenia (10). The use of atypical antipsychotics was, however, found to have a high occurrence of a potentially fatal blood disorder called agranulocytosis (10). Therefore, discovering effective, safe antipsychotics remains a high priority (4,9,10). The pathogenesis of schizophrenia was suggested to be related to dysfunction of the D1 receptor in the medial prefrontal cortex (mPFC), which is accompanied by D2 receptor hyperactivity in subcortical regions such as the ventral tegmental area (VTA) and the nucleus accumbens (NAc) (4). The D1 dysfunction was suggested to be responsible for the negative symptoms of schizophrenia, whereas the hyperactivity of the D2 receptor might lead to the positive symptoms of this disorder (4,5,11). Based on this hypothesis, antipsychotic Submitted May 8, 2006, and accepted for publication March 20, 2007. Address reprint requests to Prof. Hualiang Jiang and Weiliang Zhu, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, 555 Zu Chong Zhi Road, Zhangjiang Hi-Tech Park, Shanghai 201203, People’s Republic of China. Tel.: 86-21-50805873; Fax: 86-21-50807088; E-mail: hljiang@mail.shcnc.ac.cn. Editor: Ivet Bahar. 2007 by the Biophysical Society 0006-3495/07/09/1431/11 $2.00 doi: 10.1529/biophysj.106.088500 Biophysical Journal Volume 93 September 2007 1431–1441 1431
et drugs with dual effect as a DI receptor agonist and a D2 MODELING AND SIMULATION METHODS receptor antagonist would provide a new treatment for psy- strategy of modeling and simulation chotic diseases(4) An active compound, namely (-H-stepholidine(SPD), Experiments demonstrated that the sequences of aminergic receptors, viz. isolated from the Chinese herb Stephania, is to date the only dopamine, a-adrer drug with a dual effect as a DI receptor agonist and a D2 conservative within the TM domains, which dictate the common ligand- receptor antagonist(4). SPD has high affinity for DI- and suggest that different binding interactions of ligands to one receptor lead to D2-like receptors but low affinity for 5-HT2 receptors and different receptor conformations (10). The small cha a2-adrenoceptors(4). The dual action of SPD has been dem affect its interactions with the receptor and, hence, the receptor activa onstrated in both cortical and subcortical structures, includ tion(10). Thus, we faced two difficulties: 1), how to construct 3D models of the DI and D2 receptors, and 2), how to identify the binding conforma ing the mPFC. NAc, VTA, and basal ganglia dopamine tions of the two receptors with SPD. To solve these problems, we integrated systems(4). Moreover, clinic studies showed that SPD is several modeling and simulation methods in this study. The computational superior to perphenazine in antipsychotic efficacy. Unlike pipeline is outlined in Fig. SI in the Supplementary Material. Briefly.the perphenazine, SPd does not induce any EPSs(4), making computational flow is as follows SPD an attractive compound in studying the dual action 1. The 3D models of the DI and D2 receptors were constructed using a mechanism of a chemical for developing novel antipsychotic homology-modeling approach based on the x-ray crystal structure of drugs. However, the molecular basis of the dual action of bovine rhode SPD is obscure to date. Thus, exploring the dual action 2. Two Io-ns MD simulations were carried out respectively on the con- mechanism of SPD against DI and D2 receptors at the (palmitoyloleoylphosphatidylcholine) bilayer. atomic level is the first step toward developing more superior 3. To find the probable SPD binding conformations of DI and D2 recep- antipsychotic agents. Obviously, the 3D structures of DI and tors, SPD was docked into numerous minimum-energy conformations D2 receptors are essential to figuring out the dual action isolated from the MD trajectories of the two receptors by using the me mechanism. Unfortunately, these 3d structures are not cur- lecular docking approach, includin rently available energy between SPD and each selected conformation of the receptors. It has been widely recognized that molecular modeling Conformations of the two receptors with lowest binding free energies to SPD were selected as the initial structures for further simulations and simulation are an excellent complement to experiments 4. Two additional 10-ns MD simulations were conducted on the initial in explaining experimental results and providing clues for structures of SPD-DI and SPD-D2 complexes, resulting from the above further experiments. Furthermore, it may reveal information docking calculations, embedded in a hydrated POPC bilayer. that is not accessible by experiments(12). Recently, grea 5. The ligand-receptor interaction and the dual action mechanism of SPD were explored by analyzing all of the modeling and simulation results. success has been achieved in the field of structure predic tion of GPCRs(13-15). For instance, some algorithms ar apable of predicting the transmembrane(TM)structures (15-20). However, constructing appropriate conformational Homology modeling for the 3D models of D1 and states of the extracellular or intracellular domains of GPCRs D2 receptors is still a tough job. Molecular dynamics(MD) simulations Bovine rhodopsin, a member of the GPCR superfamily, was structur (21), taking advantage of iteratively tracking the trajectory determined at a resolution of 2.80 A(Protein Data Bank entry 1F88)(I f conformational change could be used to simulate the This x-ray structure of GPCRs provides a solid template for modeling the 3D binding process, to explore the binding conformation, and to structures of other GPCRs () To obtain a reliable map the binding mechanism at molecular and atomic levels pereceptorsdownloadedfromhttp://www. gpcr. org//tm were used. Residues were numbered according to the gener- To the best of our knowledge, no long timescale Md alized numbering scheme proposed by Ballesteros and Weinstein(22).To simulation has been carried out so far on drs, all this facilitate the comparison among the aligned residues in various GPCRs, the motivates us to carry out a computational study on DI and x 50 and residues within a given tM are then indexed relative to the 50 SPD on these two receptors. Thus, homology modeling, The MODELLER program encoded in Insightll (23)was employed to automated molecular docking, and MD simulations have assemble the 3D models of the a-helix bundles of DI and D2 receptors by been integrated in this study. A homology modeling using the x-ray crystal structure of bovine rhodopsin as a template. The approach was used to construct 3D structures of DI and FASTA program (24)was used to identify sequence homology through an D2 receptors. Automated docking was used to find possible in-house database(15)containing 700 loops and proteins with medium to high sequence identity. Clustalw (25)was then used to determine the binding sites for SPD against the receptors. The MD sim- fragments that have higher homology with the ulations, performed in a fully hydrated lipid bilayer envi- A reasonable fragment conformation was chosen from the top 10 candidates ronment, were carried out to address the following questions at have the lowest root mean-square(RMs) values and considerable How does SPD bind to the DI and D2 receptors? How does geometrical compatibility. The conserved disultide bond between residues SPD regulate the signal transe n of the DI receptor? Cys-3. 25 at the beginning of TM Ill and Cys_EL-2 in the middle of extracellular loop 2(EL-2)was also created. Energy minimizations of the What is the fund tal basis lying the dual action of models were carried out using the Discover module encoded in Insightll SPD? with the same parameters as that of our previous studies(12). Biophysical Journal 93(5)1431-1441
drugs with dual effect as a D1 receptor agonist and a D2 receptor antagonist would provide a new treatment for psychotic diseases (4). An active compound, namely ()–stepholidine (SPD), isolated from the Chinese herb Stephania, is to date the only drug with a dual effect as a D1 receptor agonist and a D2 receptor antagonist (4). SPD has high affinity for D1- and D2-like receptors but low affinity for 5-HT2 receptors and a2-adrenoceptors (4). The dual action of SPD has been demonstrated in both cortical and subcortical structures, including the mPFC, NAc, VTA, and basal ganglia dopamine systems (4). Moreover, clinic studies showed that SPD is superior to perphenazine in antipsychotic efficacy. Unlike perphenazine, SPD does not induce any EPSs (4), making SPD an attractive compound in studying the dual action mechanism of a chemical for developing novel antipsychotic drugs. However, the molecular basis of the dual action of SPD is obscure to date. Thus, exploring the dual action mechanism of SPD against D1 and D2 receptors at the atomic level is the first step toward developing more superior antipsychotic agents. Obviously, the 3D structures of D1 and D2 receptors are essential to figuring out the dual action mechanism. Unfortunately, these 3D structures are not currently available. It has been widely recognized that molecular modeling and simulation are an excellent complement to experiments in explaining experimental results and providing clues for further experiments. Furthermore, it may reveal information that is not accessible by experiments (12). Recently, great success has been achieved in the field of structure prediction of GPCRs (13–15). For instance, some algorithms are capable of predicting the transmembrane (TM) structures (15–20). However, constructing appropriate conformational states of the extracellular or intracellular domains of GPCRs is still a tough job. Molecular dynamics (MD) simulations (21), taking advantage of iteratively tracking the trajectory of conformational change, could be used to simulate the binding process, to explore the binding conformation, and to map the binding mechanism at molecular and atomic levels. To the best of our knowledge, no long timescale MD simulation has been carried out so far on DRs. All this motivates us to carry out a computational study on D1 and D2 receptors to explore the mechanism of the dual action of SPD on these two receptors. Thus, homology modeling, automated molecular docking, and MD simulations have been integrated in this study. A homology modeling approach was used to construct 3D structures of D1 and D2 receptors. Automated docking was used to find possible binding sites for SPD against the receptors. The MD simulations, performed in a fully hydrated lipid bilayer environment, were carried out to address the following questions: How does SPD bind to the D1 and D2 receptors? How does SPD regulate the signal transduction of the D1 receptor? What is the fundamental basis underlying the dual action of SPD? MODELING AND SIMULATION METHODS Strategy of modeling and simulation Experiments demonstrated that the sequences of aminergic receptors, viz. dopamine, a-adrenergic, b-adrenergic, and serotonin receptors, are highly conservative within the TM domains, which dictate the common ligandbinding sites of these receptors (22). However, experimental data also suggest that different binding interactions of ligands to one receptor lead to different receptor conformations (10). The small changes in ligand structure may affect its interactions with the receptor and, hence, the receptor activation (10). Thus, we faced two difficulties: 1), how to construct 3D models of the D1 and D2 receptors, and 2), how to identify the binding conformations of the two receptors with SPD. To solve these problems, we integrated several modeling and simulation methods in this study. The computational pipeline is outlined in Fig. S1 in the Supplementary Material. Briefly, the computational flow is as follows: 1. The 3D models of the D1 and D2 receptors were constructed using a homology-modeling approach based on the x-ray crystal structure of bovine rhodopsin. 2. Two 10-ns MD simulations were carried out respectively on the constructed D1 and D2 receptor models, embedded in a hydrated POPC (palmitoyloleoylphosphatidylcholine) bilayer. 3. To find the probable SPD binding conformations of D1 and D2 receptors, SPD was docked into numerous minimum-energy conformations isolated from the MD trajectories of the two receptors by using the molecular docking approach, including prediction of the binding free energy between SPD and each selected conformation of the receptors. Conformations of the two receptors with lowest binding free energies to SPD were selected as the initial structures for further simulations. 4. Two additional 10-ns MD simulations were conducted on the initial structures of SPD-D1 and SPD-D2 complexes, resulting from the above docking calculations, embedded in a hydrated POPC bilayer. 5. The ligand-receptor interaction and the dual action mechanism of SPD were explored by analyzing all of the modeling and simulation results. Homology modeling for the 3D models of D1 and D2 receptors Bovine rhodopsin, a member of the GPCR superfamily, was structurally determined at a resolution of 2.80 A˚ (Protein Data Bank entry 1F88) (17). This x-ray structure of GPCRs provides a solid template for modeling the 3D structures of other GPCRs (7). To obtain a reliable sequence alignment, 64 sequences of dopamine-type receptors downloaded from http://www. gpcr.org/7tm were used. Residues were numbered according to the generalized numbering scheme proposed by Ballesteros and Weinstein (22). To facilitate the comparison among the aligned residues in various GPCRs, the most conserved residue in transmembrane X is given the index number X.50, and residues within a given TM are then indexed relative to the ‘‘50’’ position. The MODELLER program encoded in InsightII (23) was employed to assemble the 3D models of the a-helix bundles of D1 and D2 receptors by using the x-ray crystal structure of bovine rhodopsin as a template. The FASTA program (24) was used to identify sequence homology through an in-house database (15) containing 700 loops and proteins with medium to high sequence identity. ClustalW (25) was then used to determine the fragments that have higher homology with the loops of D1 and D2 receptors. A reasonable fragment conformation was chosen from the top 10 candidates that have the lowest root mean-square (RMS) values and considerable geometrical compatibility. The conserved disulfide bond between residues Cys-3.25 at the beginning of TM III and Cys_EL-2 in the middle of extracellular loop 2 (EL-2) was also created. Energy minimizations of the models were carried out using the Discover module encoded in InsightII with the same parameters as that of our previous studies (12). 1432 Fu et al. Biophysical Journal 93(5) 1431–1441
D1 Agonist and D2 Antagonist Dual Effect of SPD 1433 Molecular dynamics simulations binding sites of the energy-minimum conformations of the DI and D2 receptors derived from the MD trajectories. The MD simulations were performed using the GrOMACS package ver sion3.1.4withtheGromoS96forcefield(26,27)(www.gromacs.org) The molecular topology file for SPD was generated with PRODRG (28) http://davapcl.bioch.dundeeac.uk/programs/prodrg/prodrg.html).Thepa RESULTS tial atomic charges of SPD were determined by using the CHelpG method (29) implemented in the Gaussian98 program(30)with the DFT/B3LYP/ 3D structures of D1 and D2 receptors unliganded D2, SPD-DI, and SPD-D2 complexes)were embedded in a Sequence alignment(Fig. S2 in the Supplementary Material) hydrated POPC lipid bilayer. The procedure and parameters for constructing indicates that the sequence identity and similarity are 21.8% the receptor/hydrated POPC systems are similar to those used in previous and 47.8%, respectively, for the TMs between rhodopsin membrane protein simulations(31-33). Fig. I shows, taking the SPD-D1 and the Dl receptor, 26.1% and 54.4% for the TMs between complex as an example, the structural mode of the receptor/PoPC/water rhodopsin and D2 receptor, and 44.5% and 66.4% between the four receptor/hydrated POPC systems first for all water molecules to the DI and D2 receptors themselves. Based on the high re their poor contacts with protein atoms, then for the whole system homology revealed by sequence alignments, the 3D models Im force was <10.00 kcal/mol- of the DI and D2 receptors were assembled taking the x-ray The solvent (water and PoPC) molecules of each initial system were crystal structure of bovine rhodopsin as a template.The PROCHECK (38)statistics showed that 90% of the residues molecules were constrained at 10. 50. 100. 200 and 298 K Afterward. each in both the DI and D2 models were in either the most favored system was equilibrated for 250 ps without any constraints To maintain the or in the additionally allowed regions of the ramachandran ystems at a constant temperature of 300 K, the Berendsen thermostat (34) map(Fig. S3 in the Supplementary Material), suggesting was applied using a coupling time of 0. 1 ps for the bulk water and PoPC. that the overall main chain and side chain structures are The values of the anisotropic isothermal compressibility were set to 45x all reasonable. The WHATIF (9)validation shows accept- 10-545x 10-5, 4.5 x 10-5,0,0, 0 bar"lfor xx, yy, zz, xy/yx, xz/zx, and able RMS Z-scores. All of these data suggest that the models yzizy components, respectively, for water and POPC simulations. The obtained through our homology modeling are reasonable lengths of all bonds, including those to hydrogen atoms, were constrained by Two kinds of interaction networks are observed in di and the LINCS algorithm(26). Electrostatic interactions between charged groups withing A were calculated explicitly, and long-range electrostatic interac- D2, 1.e, aromatic residue clusters(Fig. S4)and hydrogen tions were calculated using the particle mesh Ewald method (35)with grid bond (H-bond)network (Table SI). The presence of hydro- width of 1.2 A and a fourth-order spline interpolation. A cutoff distance of phobic cluster and H-bond interaction with TM Ill and TM 14 A was applied for the Lennard-Jones interactions. Numerical integratio vi in di could allow relative movement of tm vi. the of the equations of motion used a time step of 2 fs with atomic coordinates specific H-bonds cluster at the bottom of TM Ill and TM Vis saved every I ps for later analysis. To neutralize the modeled systems, 13, e added to the molecular systems of the free D1, D2, observed consisting of R-3.50 and E-6.30. It is conserved in both DI and D2 receptors, and the DRY motif of TM Ill may imulations were performed on these systems under the periodic boundary be relevant to the activation or inactivation of the receptor conditions in the nPt canonical ensemb The extracellular loop 1(EL-1)and extracellular loop 3 (EL-3)are typically short in all aminergic GPCRs. In con- trast,EL-2 is significantly longer and may reach into the active Molecular docking and binding energy calculation site crevice and form a lid over the bound ligand. There is a highly conserved disulfide bond between the conserved The geometry of the SPD ligand was built based on its crystal structure (6) Cys_e2 at the middle of EL-2 and Cy ys-3.25 at the beginning features the r configuration and a half-chair conformation of rings b and C. of TMIl. Thus in DI and D2 receptors, the stretches of only The dihedral angle between rings A and D is 1580. The molecular docking 4-5 residues between Cys_e2 and extracellular end of TM5 program AutoDock3. 05(37)was employed to probe the possible SPD. are in an extended state to reach two sides. It is noteworthy FIGURE 1 Overview of the SPD-DI/POPC/water sys- bilayer lipids are labeled H; and the lipid carbonyl oxygen atoms, choline N atoms, and P atoms are represented in a pace fill model (labeled P) and water as w: the 3-D size f the whole system was also labeled. The SPD in bindi Amplified binding SPD site was amplified and shown in a stick model. Biophysical Joumal 93(5)1431-1441
Molecular dynamics simulations The MD simulations were performed using the GROMACS package version 3.1.4 with the GROMOS96 force field. (26,27) (www.gromacs.org). The molecular topology file for SPD was generated with PRODRG (28) (http://davapc1.bioch.dundee.ac.uk/programs/prodrg/prodrg.html). The partial atomic charges of SPD were determined by using the CHelpG method (29) implemented in the Gaussian98 program (30) with the DFT/B3LYP/ 6-311G** basis set. For MD simulations, the four models (unliganded D1, unliganded D2, SPD-D1, and SPD-D2 complexes) were embedded in a hydrated POPC lipid bilayer. The procedure and parameters for constructing the receptor/hydrated POPC systems are similar to those used in previous membrane protein simulations (31–33). Fig. 1 shows, taking the SPD-D1 complex as an example, the structural model of the receptor/POPC/water systems. Before MD simulations, energy minimizations were performed on the four receptor/hydrated POPC systems first for all water molecules to remove their poor contacts with protein atoms, then for the whole system until the maximum force was ,10.00 kcal/mol-A˚ . The solvent (water and POPC) molecules of each initial system were equilibrated with protein structures by constraining the solute (D1 or D2) at 300 K for 20 ps. Then the protein was equilibrated for 5 ps and the solvent molecules were constrained at 10, 50, 100, 200, and 298 K. Afterward, each system was equilibrated for 250 ps without any constraints. To maintain the systems at a constant temperature of 300 K, the Berendsen thermostat (34) was applied using a coupling time of 0.1 ps for the bulk water and POPC. The pressure was maintained by coupling to a reference pressure of 1 bar. The values of the anisotropic isothermal compressibility were set to 4.5 3 105 , 4.53 105 , 4.5 3 105 , 0, 0, 0 bar1 for xx, yy, zz, xy/yx, xz/zx, and yz/zy components, respectively, for water and POPC simulations. The lengths of all bonds, including those to hydrogen atoms, were constrained by the LINCS algorithm (26). Electrostatic interactions between charged groups within 9 A˚ were calculated explicitly, and long-range electrostatic interactions were calculated using the particle mesh Ewald method (35) with a grid width of 1.2 A˚ and a fourth-order spline interpolation. A cutoff distance of 14 A˚ was applied for the Lennard-Jones interactions. Numerical integration of the equations of motion used a time step of 2 fs with atomic coordinates saved every 1 ps for later analysis. To neutralize the modeled systems, 13, 9, 12, and 8 Cl ions were added to the molecular systems of the free D1, D2, SPD-D1, and SPD-D2 complexes, respectively. Finally, four 10-ns MD simulations were performed on these systems under the periodic boundary conditions in the NPT canonical ensemble. Molecular docking and binding energy calculation The geometry of the SPD ligand was built based on its crystal structure (36) and optimized at the DFT/B3LYP/6-311G** level. This protonated structure features the R configuration and a half-chair conformation of rings B and C. The dihedral angle between rings A and D is 158. The molecular docking program AutoDock3.05 (37) was employed to probe the possible SPDbinding sites of the energy-minimum conformations of the D1 and D2 receptors derived from the MD trajectories. RESULTS 3D structures of D1 and D2 receptors Sequence alignment (Fig. S2 in the Supplementary Material) indicates that the sequence identity and similarity are 21.8% and 47.8%, respectively, for the TMs between rhodopsin and the D1 receptor, 26.1% and 54.4% for the TMs between rhodopsin and D2 receptor, and 44.5% and 66.4% between the D1 and D2 receptors themselves. Based on the high homology revealed by sequence alignments, the 3D models of the D1 and D2 receptors were assembled taking the x-ray crystal structure of bovine rhodopsin as a template. The PROCHECK (38) statistics showed that 90% of the residues in both the D1 and D2 models were in either the most favored or in the additionally allowed regions of the Ramachandran map (Fig. S3 in the Supplementary Material), suggesting that the overall main chain and side chain structures are all reasonable. The WHATIF (39) validation shows acceptable RMS Z-scores. All of these data suggest that the models obtained through our homology modeling are reasonable. Two kinds of interaction networks are observed in D1 and D2, i.e., aromatic residue clusters (Fig. S4) and hydrogen bond (H-bond) network (Table S1). The presence of hydrophobic cluster and H-bond interaction with TM III and TM VII in D1 could allow relative movement of TM VI. The specific H-bonds cluster at the bottom of TM III and TM V is observed consisting of R-3.50 and E-6.30. It is conserved in both D1 and D2 receptors, and the DRY motif of TM III may be relevant to the activation or inactivation of the receptor. The extracellular loop 1 (EL-1) and extracellular loop 3 (EL-3) are typically short in all aminergic GPCRs. In contrast, EL-2 is significantly longer and may reach into the active site crevice and form a lid over the bound ligand. There is a highly conserved disulfide bond between the conserved Cys_e2 at the middle of EL-2 and Cys-3.25 at the beginning of TM III. Thus in D1 and D2 receptors, the stretches of only 4–5 residues between Cys_e2 and extracellular end of TM5 are in an extended state to reach two sides. It is noteworthy FIGURE 1 Overview of the SPD-D1/POPC/water system. D1 is shown in ribbon; the hydrophobic chains of bilayer lipids are labeled H; and the lipid carbonyl oxygen atoms, choline N atoms, and P atoms are represented in a space fill model (labeled P), and water as W; the 3-D size of the whole system was also labeled. The SPD in binding site was amplified and shown in a stick model. D1 Agonist and D2 Antagonist Dual Effect of SPD 1433 Biophysical Journal 93(5) 1431–1441
Fu et al that there are 13 residues in the EL-2 loop of the D2 receptor, whole DI receptor and peak-minimum distance changes of whereas 26 residues exist in the EL-2 loop of the DI K-167-E-302(Fig 4 A)occurred three times during our 10- eptor. The main difference in EL-2 mainly comes from ns MD simulations. Quite similar open-closed conforma the upstream preceding conserved Cys_e2(Fig. S2), which is tional changes were also seen for the free D2(Fig. 3 B) far away from the active site crevice. In addition, there These dynamic properties suggest that there are at least two six residues longer in the EL-3 of the DI receptor than that of distinct conformations for both DI or D2, and the sponta the D2 receptor. Though the lengths of EL-2 and EL-3 are neous oscillation between these two conformational states different, they are flexible and function as the conformational appears to be an intrinsic phenomenon. Our observation is in switch to catch the outcoming ligands(see the discussion in agreement with experimental mutagenesis studies in which it the next paragraph) was postulated that GPCRs exist in an equilibrium between two interchangeable conformational states(40, 41). One of the obvious advantages of this type of dynamical behavior of Open-closed conformational changes of free the binding site entrances(open-closed) is that it enables the D1 and D2 receptors to easily capture different ligands(agonists and antagonists)or different conformations of a ligand, as will be The MD simulations of the free DI and D2 systems showed discussed later that, for all of the systems, the temperature, mass density, Tracing the MD trajectory of both the DI and D2 receptors and volume are relatively stable after 2 ns. About then, the revealed a portion of the atomic-level mechanism of the fuctuation scale became much smaller for both the RMs open-closed motion. A hydrogen-bond(H-bond) network deviations of the Ca atoms and potential energies(Fig. 2)of exists between the charged side chains of K-167 and E-302 the two simulation systems, indicating that the molecular At the beginning of MD simulation, the positively charged ystems were well behaved thereafter. side chain of K-167 forms direct H-bonds and two indirect When analyzing the MD trajectories, a striking"open-(water-bridged) H-bonds with the negatively charged side closed"conformational change is observed from the extracel- chain of E-302(Fig. 5 A). As the simulation proceeds,the lular side for both the DI and D2 receptors, as demonstrated side chain of the E-302 turns away from K-167, breaking the Fig 3. Further examination reveals that such an"open- H-bond network. At-2630 ps, all of the H-bonds between closed'event can be attributed mainly to large conforma- K-167 and E-302 are gone. Meanwhile, K-167 gradually tional changes in the EL-2 and EL-3, coupled with the approaches D-173 of EL-2, forming a new H-bond network connecting TM helices. Typically, EL-2 and EL-3 act as via water bridges(Fig. 5 B)which lasts -I ns. K-167-E-302 lips'"of mouth-like entrances of the binding sites of these switches back to the pairing state similar to that in Fig. 5 A two receptors. This can be demonstrated by the fluctuations and the H-bond network is restored. K-167 changes partners of critical distance between K-167 at EL-2 and E-302 at EL- frequently, either E-302 or D-173, staying with D-173 when 3(K-167-E-302 pairing), which is an objective monitor of the"mouth"is open and with E-302 when the"mouth"is he width of the entrance. The open and closed conforma- closed. Therefore, the H-bond network between K-167 and tions are shown in the dotted lines in Fig. 2. Taking DI as an E-302 seems to act as a"bolt"its position control example, the opening event at -2630 ps for the"mouth whether the"mouth"is open or closed orresponds to changes in the K-167-E-302 pairing distance (Fig 4 A), which reaches one of the "peak" points as shown in Fig. 6. As time proceeds, the"mouth"becomes gradually Binding conformations of SPD in the closed at -4210 ps and the distance between K-167-E-302 D1 and D2 receptors decreases to a local minimum. Such synchronization To identify the most probable binding conformations of SPD between the"open-closed"conformational changes of the to the DI and D2 receptors, 10 potential binding conformation 457000 FIGURE 2 Energy changes along MD mulations. (A)DI simulation system. (B)D2 simulation system. Ten conforma- 45900 Most probably active ons shown in dotted line were identified 243m0 conformation of 2 460000 nd obvious geometrical difference among 346100 different conformations. Among the 4453000 the conformations at 2630 ps in the MD 463000 ajectory of the unliganded DI and that at 310 ps in the MD trajectory of ur 020040060080001000 0 2000 4000 6000 8000 10000 liganded D2 have the lowest binding time(ps) energies toward SPD and were selected as the conformations most probably active. Biophysical Journal 93(5)1431-1441
that there are 13 residues in the EL-2 loop of the D2 receptor, whereas 26 residues exist in the EL-2 loop of the D1 receptor. The main difference in EL-2 mainly comes from the upstream preceding conserved Cys_e2 (Fig. S2), which is far away from the active site crevice. In addition, there are six residues longer in the EL-3 of the D1 receptor than that of the D2 receptor. Though the lengths of EL-2 and EL-3 are different, they are flexible and function as the conformational switch to catch the outcoming ligands (see the discussion in the next paragraph). Open-closed conformational changes of free D1 and D2 The MD simulations of the free D1 and D2 systems showed that, for all of the systems, the temperature, mass density, and volume are relatively stable after 2 ns. About then, the fluctuation scale became much smaller for both the RMS deviations of the Ca atoms and potential energies (Fig. 2) of the two simulation systems, indicating that the molecular systems were well behaved thereafter. When analyzing the MD trajectories, a striking ‘‘openclosed’’ conformational change is observed from the extracellular side for both the D1 and D2 receptors, as demonstrated in Fig. 3. Further examination reveals that such an ‘‘openclosed’’ event can be attributed mainly to large conformational changes in the EL-2 and EL-3, coupled with the connecting TM helices. Typically, EL-2 and EL-3 act as ‘‘lips’’ of mouth-like entrances of the binding sites of these two receptors. This can be demonstrated by the fluctuations of critical distance between K-167 at EL-2 and E-302 at EL- 3 (K-167-E-302 pairing), which is an objective monitor of the width of the entrance. The open and closed conformations are shown in the dotted lines in Fig. 2. Taking D1 as an example, the opening event at ;2630 ps for the ‘‘mouth’’ corresponds to changes in the K-167-E-302 pairing distance (Fig. 4 A), which reaches one of the ‘‘peak’’ points as shown in Fig. 6. As time proceeds, the ‘‘mouth’’ becomes gradually closed at ;4210 ps and the distance between K-167-E-302 decreases to a local minimum. Such synchronization between the ‘‘open-closed’’ conformational changes of the whole D1 receptor and peak-minimum distance changes of K-167-E-302 (Fig. 4 A) occurred three times during our 10- ns MD simulations. Quite similar open-closed conformational changes were also seen for the free D2 (Fig. 3 B). These dynamic properties suggest that there are at least two distinct conformations for both D1 or D2, and the spontaneous oscillation between these two conformational states appears to be an intrinsic phenomenon. Our observation is in agreement with experimental mutagenesis studies in which it was postulated that GPCRs exist in an equilibrium between two interchangeable conformational states (40,41). One of the obvious advantages of this type of dynamical behavior of the binding site entrances (open-closed) is that it enables the receptors to easily capture different ligands (agonists and antagonists) or different conformations of a ligand, as will be discussed later. Tracing the MD trajectory of both the D1 and D2 receptors revealed a portion of the atomic-level mechanism of the open-closed motion. A hydrogen-bond (H-bond) network exists between the charged side chains of K-167 and E-302. At the beginning of MD simulation, the positively charged side chain of K-167 forms direct H-bonds and two indirect (water-bridged) H-bonds with the negatively charged side chain of E-302 (Fig. 5 A). As the simulation proceeds, the side chain of the E-302 turns away from K-167, breaking the H-bond network. At ;2630 ps, all of the H-bonds between K-167 and E-302 are gone. Meanwhile, K-167 gradually approaches D-173 of EL-2, forming a new H-bond network via water bridges (Fig. 5 B) which lasts ;1 ns. K-167-E-302 switches back to the pairing state similar to that in Fig. 5 A, and the H-bond network is restored. K-167 changes partners frequently, either E-302 or D-173, staying with D-173 when the ‘‘mouth’’ is open and with E-302 when the ‘‘mouth’’ is closed. Therefore, the H-bond network between K-167 and E-302 seems to act as a ‘‘bolt’’—its position controls whether the ‘‘mouth’’ is open or closed. Binding conformations of SPD in the D1 and D2 receptors To identify the most probable binding conformations of SPD to the D1 and D2 receptors, 10 potential binding conformation FIGURE 2 Energy changes along MD simulations. (A) D1 simulation system. (B) D2 simulation system. Ten conformations shown in dotted line were identified with two criteria: lower potential energy and obvious geometricaldifferenceamong different conformations. Among them, the conformations at 2630 ps in the MD trajectory of the unliganded D1 and that at 3310 ps in the MD trajectory of unliganded D2 have the lowest binding energies toward SPD and were selected as the conformations most probably active. 1434 Fu et al. Biophysical Journal 93(5) 1431–1441
D1 Agonist and D2 Antagonist Dual Effect of SPD FIGURE ns of the DI (A)and D2 (B)receptor during the open-closed process of the inding site as observed from MD simu- lations. a dashed arrow line indicates the position of the"mouth, and a double dashed line shows the*"mouth 緲嚕輬 pening. All of the conformations are hown in the style of molecular surfa plored by electrostatic potential. The ed color stands for negative electrostatic candidates of DRs were identified from each of the MD Dual action mechanism of SPD with dopamine rajectories of the unliganded DI and D2 receptors(dotted D1 and D2 receptors lines in Fig. 2)with two criteria: lower potential energy and The structural determinant of pharmacological specificity obvious geometrical difference among different conforma- of SPD ons. Then spd was docked into the binding cavities of these 10 conformations of each receptor and the binding The dynamics of the interactions between SPD and DI or d2 energies were estimated by means of Auto Dock3.05(37). was investigated to understand the agonistic and antagonistic Both the binding mode and binding affinity between SPD mechanism of SPD ind DRs were used to guide the selection of the conforma The distances between the key residues in the active sites ion that is most likely to be active. From among 10 were monitored and used to delineate the dynamical change candidates, the conformations were selected as the possible of the active site among the bound complexes. Most dis- conformation candidates if their main binding modes are in tances fluctuate a little bit in complexes SPD-DI. It shows agreement with known mutagenesis experiments. From these that SPD packs well with the side chains of residues in the possible conformation candidates, the conformation having binding cavity of DI. In the SPD-D2 complex, there is an the lowest binding energy toward SPD was selected as the interesting conformational change around the ligand. The conformations most likely to be active for each DR. The side chains of F-6.52 and H-6.55 come near the aromatic predicted lowest binding free energy of SPD for Dl is.8 rings A and D of SPD, apparently induced by hydrophobic kcal/mol and.7 kcal/mol, respectively. Although there interactions among these groups and the bending of rings A are deviations between the experimental value(-10.8 kcal/ and d toward side chains of F-6.52 and H-6.55. There is also mol and.6 kcal/mol)(42-45)and predicted data, the an electrostatic attraction between the protonated side chain general trend observed in predicated binding free energies is of H-6.55 and the electron-rich group of ring A and atom that SPd binds to DI more strongly than to D2. The errors in O2 of SPD(Fig. 6 D). The conformational bending of SPD predicting binding affinities mainly come from the homol- and the aromatic packing with the side chains of F-6.52 ogy models of DRs and the imperfect empirical parameters and H-6.55 appear to hold SPD against helix VI, constrain- in AutoDock 3.05 for estimating binding affinity ing it after SPD binding. In comparison, there is no such E Away FIGURE 4 Distance fluctuations during the MD simulations. (A)Distance from Na B他 of K-167(EL-2)to Cs of E-302(EL-3)of the DI receptor; (B)distance from Cy of 08 D-178(EL-2)to the centroid of the are 6H matic side chain of Y-408(EL-3)of the D2 04 receptor. The open and closed conforma- state open state tions in Fig 3 are shown in the dotted lines me("0 800010000 200040006000080001000inhg.2. B Biophysical Joumal 93(5)1431-1441
candidates of DRs were identified from each of the MD trajectories of the unliganded D1 and D2 receptors (dotted lines in Fig. 2) with two criteria: lower potential energy and obvious geometrical difference among different conformations. Then SPD was docked into the binding cavities of these 10 conformations of each receptor and the binding energies were estimated by means of AutoDock3.05 (37). Both the binding mode and binding affinity between SPD and DRs were used to guide the selection of the conformation that is most likely to be active. From among 10 candidates, the conformations were selected as the possible conformation candidates if their main binding modes are in agreement with known mutagenesis experiments. From these possible conformation candidates, the conformation having the lowest binding energy toward SPD was selected as the conformations most likely to be active for each DR. The predicted lowest binding free energy of SPD for D1 is 12.8 kcal/mol and 11.7 kcal/mol, respectively. Although there are deviations between the experimental value (10.8 kcal/ mol and 9.6 kcal/mol) (42–45) and predicted data, the general trend observed in predicated binding free energies is that SPD binds to D1 more strongly than to D2. The errors in predicting binding affinities mainly come from the homology models of DRs and the imperfect empirical parameters in AutoDock 3.05 for estimating binding affinity. Dual action mechanism of SPD with dopamine D1 and D2 receptors The structural determinant of pharmacological specificity of SPD The dynamics of the interactions between SPD and D1 or D2 was investigated to understand the agonistic and antagonistic mechanism of SPD. The distances between the key residues in the active sites were monitored and used to delineate the dynamical change of the active site among the bound complexes. Most distances fluctuate a little bit in complexes SPD-D1. It shows that SPD packs well with the side chains of residues in the binding cavity of D1. In the SPD-D2 complex, there is an interesting conformational change around the ligand. The side chains of F-6.52 and H-6.55 come near the aromatic rings A and D of SPD, apparently induced by hydrophobic interactions among these groups and the bending of rings A and D toward side chains of F-6.52 and H-6.55. There is also an electrostatic attraction between the protonated side chain of H-6.55 and the electron-rich group of ring A and atom O2 of SPD (Fig. 6 D). The conformational bending of SPD and the aromatic packing with the side chains of F-6.52 and H-6.55 appear to hold SPD against helix VI, constraining it after SPD binding. In comparison, there is no such FIGURE 3 Representative conformations of the D1 (A) and D2 (B) receptors during the open-closed process of the binding site as observed from MD simulations. A dashed arrow line indicates the position of the ‘‘mouth’’, and a double arrow dashed line shows the ‘‘mouth’’ opening. All of the conformations are shown in the style of molecular surface colored by electrostatic potential. The red color stands for negative electrostatic potential, and the blue for positive potential. FIGURE 4 Distance fluctuations during the MD simulations. (A) Distance from Nz of K-167 (EL-2) to Cd of E-302 (EL-3) of the D1 receptor; (B) distance from Cg of D-178 (EL-2) to the centroid of the aromatic side chain of Y-408 (EL-3) of the D2 receptor. The open and closed conformations in Fig. 3 are shown in the dotted lines in Fig. 2. D1 Agonist and D2 Antagonist Dual Effect of SPD 1435 Biophysical Journal 93(5) 1431–1441