PHYSICAL CHEMISTRY B一 pubs. acs. org/JPCB Molecular Insights into the D1R Agonist and d2R/D3R Antagonist Effects of the Natural Product (-) -Stepholidine: Molecular Modeling and dynamics simulations Bian.ts Wei Li. t, peng Du, Kun Qian Yu, *+ and Wei Fu*+,? Department of Medicinal Chemistry Key Laboratory of Smart Drug Delivery, Ministry of Education PLA, School of Pharmacy, Fudan University, Shanghai 201203, China Drug Discovery and Design Center, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China S Supporting Information ABSTRACT:(-)-Stepholidine(1SPD), an active ingredier of the Chinese herb Stephania, is the first compound found to ave a dual function as a dopamine receptor DI agonist and ary dynamical beh of DI and d2R and their interaction modes with -spd were a investigated in our previous study. Recently, the pharmaco- logical effect of 1-SPD on D3R was elucidated as an antagonist. This new discovery in combination with the explosion of structural biology in GPCR superfamily prompted us to I-SPD perform a more comprehensive investigation on the special pharmacological profiles of 1-SPD on dopamine receptors. In this study, the integration of homology modeling, automated molecular docking, and MD simulations was used to probe the agonistic and antagonistic mechanism of 1-SPD on DIR, D2R, and D3R. Our analyses showed that hydrogen bonding of the hydroxyl group on the d ring of I-SPD with side chain of N6.55 which, in combination with hydrophobic stacking between 13.40, F6.44 and w6.48, is the key feature to mediate the agonist effect of 1-SPD on DiR, whereas the absence of hydrophobic stacking between 13.40, F6.44, and W6. 48 in D2R and D3R excludes receptor activation. Finally, the agonistic and antagonistic mechanisms of 1-SPD and an activation model of dir were proposed on the basis of these findings. The present study could guide future experimental works on these receptors and has the significance to the design of functionally selective drugs targeting dopamine receptors. ■| NTRODUCTION these inactive-state receptors as template is also problematic. G-protein-coupled receptors( GPCRs)are the largest family of The newly available active-state structures of BAR, P1 AR, and integral membrane proteins that mediate most of physiological AzA receptor thus provide a solid basis for the elucidation of the and environmental conformational changes associated with agonist binding esponses to hormones, neurotransmitters, and environmental stimulants. GPCRs constitute the largest class of therapeutic Dopamine(DA), the endogenous ligand of dopaminergic targets as evidenced by the fact that at least one-fourth of drugs neurotransmission systems, has been associated with many on the market exert their therapeutic activity by modulating physiological functions such as fine movement coordination, hodopsin-like GPCRs. In the past 3 years, there were cognition, and emotion. DA exerts its effects by activating five remarkable advances in the crystallography of G-protein- distinct dopamine receptors(Drs) which belong to the GPCR oupled receptors (GPCRs). Highlights have included the superfamily and are classified into two subfamilies:D1-like characterization of the crystal structures of antagonist-bound (DIR and DSR)and D2-like(D2r, D3R, and D4R) based on and agonist-bound GPCRs: the human B2 adrenergic receptor pharmacological and functional characteristics. It has bee (2AR),- the turkey B,AR,the human A2A adenosine established that DRs are primary targets of antipsychotic drugs receptor,",the human chemokine receptor CXCR4, and the used to treat psychomotor diseases such as schizophre human dopamine D3 receptor. The availability of such high debilitating mental illness which affects 0.5-1.5% of the resolution crystal structures can greatly help to guide the worldwide population. The pathogenesis of schizophrenia is structure-based GPCR drug discovery. However, with the suggested to be related to dysfunction of the DiR in the medial inactive-state structures of B,AR, PAR, and axa receptor medicinal chemists could probably only be able to develop Received: May ligands that stabilize the inactive conformation. The agonisti Revised: 122 nechanism of GPCRs based on homology modeling taking Published: June 15, 2012 ACS Publications 2012 American Chemical Society 121 dx. dolora/10.1021/p30492351Phys. Chem. 62012116,8121-8130
Molecular Insights into the D1R Agonist and D2R/D3R Antagonist Effects of the Natural Product (−)-Stepholidine: Molecular Modeling and Dynamics Simulations Bian Li,†,§ Wei Li,†,§ Peng Du,†,§ Kun Qian Yu,*,‡ and Wei Fu*,† † Department of Medicinal Chemistry & Key Laboratory of Smart Drug Delivery, Ministry of Education & PLA, School of Pharmacy, Fudan University, Shanghai 201203, China ‡ Drug Discovery and Design Center, State Key Laboratory of Drug Research, Shanghai Institute of Materia Medica, Chinese Academy of Sciences, Shanghai 201203, China *S Supporting Information ABSTRACT: (−)-Stepholidine (l-SPD), an active ingredient of the Chinese herb Stephania, is the first compound found to have a dual function as a dopamine receptor D1 agonist and D2 antagonist. The preliminary dynamical behaviors of D1R and D2R and their interaction modes with l-SPD were investigated in our previous study. Recently, the pharmacological effect of l-SPD on D3R was elucidated as an antagonist. This new discovery in combination with the explosion of structural biology in GPCR superfamily prompted us to perform a more comprehensive investigation on the special pharmacological profiles of l-SPD on dopamine receptors. In this study, the integration of homology modeling, automated molecular docking, and MD simulations was used to probe the agonistic and antagonistic mechanism of l-SPD on D1R, D2R, and D3R. Our analyses showed that hydrogen bonding of the hydroxyl group on the D ring of l-SPD with side chain of N6.55 which, in combination with hydrophobic stacking between I3.40, F6.44 and W6.48, is the key feature to mediate the agonist effect of l-SPD on D1R, whereas the absence of hydrophobic stacking between I3.40, F6.44, and W6.48 in D2R and D3R excludes receptor activation. Finally, the agonistic and antagonistic mechanisms of l-SPD and an activation model of D1R were proposed on the basis of these findings. The present study could guide future experimental works on these receptors and has the significance to the design of functionally selective drugs targeting dopamine receptors. ■ INTRODUCTION G-protein-coupled receptors (GPCRs) are the largest family of integral membrane proteins that mediate most of physiological responses to hormones, neurotransmitters, and environmental stimulants. GPCRs constitute the largest class of therapeutic targets as evidenced by the fact that at least one-fourth of drugs on the market exert their therapeutic activity by modulating rhodopsin-like GPCRs.1 In the past 3 years, there were remarkable advances in the crystallography of G-proteincoupled receptors (GPCRs). Highlights have included the characterization of the crystal structures of antagonist-bound and agonist-bound GPCRs: the human β2 adrenergic receptor (β2AR),2−5 the turkey β1AR,6,7 the human A2A adenosine receptor,8,9 the human chemokine receptor CXCR4,10 and the human dopamine D3 receptor.11 The availability of such highresolution crystal structures can greatly help to guide the structure-based GPCR drug discovery. However, with the inactive-state structures of β2AR, β1AR, and A2A receptor, medicinal chemists could probably only be able to develop ligands that stabilize the inactive conformation. The agonistic mechanism of GPCRs based on homology modeling taking these inactive-state receptors as template is also problematic. The newly available active-state structures of β2AR, β1AR, and A2A receptor thus provide a solid basis for the elucidation of the conformational changes associated with agonist binding. Dopamine (DA), the endogenous ligand of dopaminergic neurotransmission systems, has been associated with many physiological functions such as fine movement coordination, cognition, and emotion. DA exerts its effects by activating five distinct dopamine receptors (DRs) which belong to the GPCR superfamily and are classified into two subfamilies: D1-like (D1R and D5R) and D2-like (D2R, D3R, and D4R) based on their pharmacological and functional characteristics. It has been established that DRs are primary targets of antipsychotic drugs used to treat psychomotor diseases such as schizophrenia, a debilitating mental illness which affects 0.5−1.5% of the worldwide population.12 The pathogenesis of schizophrenia is suggested to be related to dysfunction of the D1R in the medial Received: May 21, 2012 Revised: June 14, 2012 Published: June 15, 2012 Article pubs.acs.org/JPCB © 2012 American Chemical Society 8121 dx.doi.org/10.1021/jp3049235 | J. Phys. Chem. B 2012, 116, 8121−8130
The Journal of Physical Chemistry B Article prefrontal cortex (mPFC), which is accompanied by D2R MATERIALS AND METHODS hyperactivity in subcortical regions such as the ventral tegmental area(VTA)and the nucleus accumbens(Nac) Homology Modeling of D1R and D2R.The quences of DiR and D2R were collected from the The dir dysfunction is believed to be responsible for the UniProtKB database (accession code: P21728 for DIR and negative symptoms of schizophrenia, whereas the D2R P14416 for D2R), and sequence similarity searches for DIR hyperactivity might lead to the positive symptoms of the and D2R were performed using the NCBI BLAST server. 26 disorder.3-15 Thus, an effective antipsychotic drug should The newly disclosed active-state structure of P2AR (PDB code simultaneously possess DiR agonistic and D2R antagonistic 3POG) was selected as the template to construct the agonistic activities . I5 conformation of DiR, while for D2R homology modeling, the ()-Stepholidine (I-SPD), an active natural compoun inactive-state structure of D3R(PDB code: 3PBL) was used isolated from the Chinese herb Stephania, is to date the only as the template as it shares a sequence identity of -78% with drug with a dual effect as a DiR agonist and D2R antagonist. [- D2R in transmembrane(Tm)regions. Sequence alignments of spd displays high affinity to DI-and D2-like receptors but low DIR and D2R with their templates were carried out us affinity to 5-HT, receptors and rare affinity to several other ClustalX 2.0.11 program. The homology modeling and loop neurotransmitter receptors. 5-17 Moreover, clinical investiga refinement were optimized with Modeler 9v4. Fifty homology tions showed that I-SPD is superior to perphenazine in models were generated for both DiR and D2R after loop antipsychotic efficacy. Unlike perphenazine, I-SPD does not refinement, and the one with the lowest dope score adopted for subsequent energy minimization. The canonical induce any enound in studying the dual action mechanism of disulfide bridge between residues C3.25(in Ballesteros- a chemical for developing novel antipsychotics. 18 The recent Weinstein numbering, a single most conserved residue among discovery that I-SPD has an antagonistic effect on D3R the class A GPCRs is designated x 50, where x is the endered both 1-SPD and D3R much more intriguing. However transmembrane helix number. all other residues on that helix are numbered relative to this conserved position)at the tip of SpD still remains obscure despite it was partly explained by our TM3, and cysteine of ECL2 in the midle of extracellular loop preliminary computational study. A more comprehensive understanding of the triple action mechanism of I-SPD against 1000 cycles of conjugate grades ycles of steepest descent and nodels were subjected to 1000 nt energy minimization. The DIR,D2R, and D3R at the atomic level is fundamental for procHEck program which analyzes the residue-by-residue developing superior antipsychotics. Obviously, the 3D geometry and overall structure quality was employed to structures of dir and D2R dispensable for computa evaluate the stereochemical quality of the energy-minimized tiona ing out the triple action mechanism. Unfortu- models of DIR and D2R. The starting structure of D3R in the nately, only the 3D structure of D3R is currently available simulation study comes from the crystal structure(PDB code It has been widely recognized that molecular modeling and 3PBL) simulation are an excellent complement to experiments Parameterization of /SPD and Molecular Docking explaining experimental results and guiding further experiment The geometry of I-SPD was built based on its crystal structure' design. Furthermore, additional information that is not using the SYBYL 6.9 program and optimized at the DFT/ accessible by experiments can be revealed by computational B3LYP/6-31G* level. It features the same configuration and approaches conformation as our previous study 20 us study 0 The partial atomic computational homology modeling methods based on high charges of I-SPD was determined with the DFT/B3LYP/6- sequence identity can be employed to predict the 3D models 31G** basis set by using the Gaussian 09 program with enough accuracy, which can subsequently be used to aid in Parameterization and topology generation of l-SPD for the understanding of protein functional mechanisms. 23 CHARMM force field were performed by using the Paratool Molecular dynamics(MD)simulations, taking advantage of plug-in of VMD 33 The automated molecular docking program iteratively tracking the trajectory of conformational change, can AutoDock4.0 with AutoDock Tools was employed to probe the ide information about the conformational properties of possible complex of 1-SPD bound with energy-minimized nolecular system and the way in which the conformation dopamine receptors(DIR, D2R, and D3R). Binding energy changes with time, such as the binding process of a ligand into a prediction was carried out using X-Score. All available protein's binding pocket, the binding conformation that the experimental data in combination with the lowest binding ligand and protein adopt, and other time-dependent properties. free energy were used to select the initial structure with the recent insightful information about the activation subsequent MD simulations mechanism of GPCR, 4 as well as the advent of the active- System Setup. The orientation of DiR, D2R, and D3R state structure of B2AR+s and the inactive-state structure of was aligned to the z-axis using Orient plug-in of VMD.Helmut D3R, we are motivated to perform a more comprehensive Grubmuller's SoLvATE program was then used to solvate the omputational study on DIR, D2R, and D3R to explore proteins followed by elimination of water molecules located underlying agonistic and antagonistic mechanisms of 1-SPD the hydrophobic protein-membrane interface. A pre-equili- which is essential for developing superior antipsychotics. Thus, brated rectangular patch of POPC membrane bilayer with homology modeling, automated molecular docking, and MD solvated lipid headgroups was generated by using the mulations are integrated in this study with the eventual aim to Membrane Builder plug-in of VMD. The partially solvated develop novel superior antipsychotic agents with less side proteins were immersed in the POPC membrane bilayer by ffects and tailor these agents to possess the designe matching the center of mass of proteins with that of membrane, properties for individual therapy espectively. The resulting placement of each receptor in 8122 dx. dolora/o.021/p30492351 Phys. Chem. B2012116,8121-813
prefrontal cortex (mPFC), which is accompanied by D2R hyperactivity in subcortical regions such as the ventral tegmental area (VTA) and the nucleus accumbens (NAc). The D1R dysfunction is believed to be responsible for the negative symptoms of schizophrenia, whereas the D2R hyperactivity might lead to the positive symptoms of the disorder.13−15 Thus, an effective antipsychotic drug should simultaneously possess D1R agonistic and D2R antagonistic activities.15 (−)-Stepholidine (l-SPD), an active natural compound isolated from the Chinese herb Stephania, is to date the only drug with a dual effect as a D1R agonist and D2R antagonist. lSPD displays high affinity to D1- and D2-like receptors but low affinity to 5-HT2 receptors and rare affinity to several other neurotransmitter receptors.15−17 Moreover, clinical investigations showed that l-SPD is superior to perphenazine in antipsychotic efficacy.18 Unlike perphenazine, l-SPD does not induce any extrapyramidal symptoms (EPS),19 making it an attractive compound in studying the dual action mechanism of a chemical for developing novel antipsychotics.18 The recent discovery that l-SPD has an antagonistic effect on D3R rendered both l-SPD and D3R much more intriguing. However, the complete molecular basis of the dual action mechanism of lSPD still remains obscure despite it was partly explained by our preliminary computational study.20 A more comprehensive understanding of the triple action mechanism of l-SPD against D1R, D2R, and D3R at the atomic level is fundamental for developing superior antipsychotics. Obviously, the 3D structures of D1R and D2R are indispensable for computationally figuring out the triple action mechanism. Unfortunately, only the 3D structure of D3R is currently available. It has been widely recognized that molecular modeling and simulation are an excellent complement to experiments in explaining experimental results and guiding further experiment design.20,21 Furthermore, additional information that is not accessible by experiments can be revealed by computational approaches.22 In the absence of experimental structures, computational homology modeling methods based on high sequence identity can be employed to predict the 3D models with enough accuracy, which can subsequently be used to aid in the understanding of protein functional mechanisms.23 Molecular dynamics (MD) simulations, taking advantage of iteratively tracking the trajectory of conformational change, can provide information about the conformational properties of molecular system and the way in which the conformation changes with time, such as the binding process of a ligand into a protein’s binding pocket, the binding conformation that the ligand and protein adopt, and other time-dependent properties. With the recent insightful information about the activation mechanism of GPCR,24,25 as well as the advent of the activestate structure of β2AR4,5 and the inactive-state structure of D3R,11 we are motivated to perform a more comprehensive computational study on D1R, D2R, and D3R to explore the underlying agonistic and antagonistic mechanisms of l-SPD, which is essential for developing superior antipsychotics. Thus, homology modeling, automated molecular docking, and MD simulations are integrated in this study with the eventual aim to develop novel superior antipsychotic agents with less side effects and tailor these agents to possess the designed properties for individual therapy. ■ MATERIALS AND METHODS Homology Modeling of D1R and D2R. The amino acid sequences of D1R and D2R were collected from the UniProtKB database (accession code: P21728 for D1R and P14416 for D2R), and sequence similarity searches for D1R and D2R were performed using the NCBI BLAST server.26 The newly disclosed active-state structure of β2AR (PDB code: 3P0G)4 was selected as the template to construct the agonistic conformation of D1R, while for D2R homology modeling, the inactive-state structure of D3R (PDB code: 3PBL)11 was used as the template as it shares a sequence identity of ∼78% with D2R in transmembrane (TM) regions. Sequence alignments of D1R and D2R with their templates were carried out using ClustalX 2.0.11 program.27 The homology modeling and loop refinement were optimized with Modeler 9v4.28 Fifty homology models were generated for both D1R and D2R after loop refinement, and the one with the lowest DOPE score28 was adopted for subsequent energy minimization. The canonical disulfide bridge between residues C3.25 (in Ballesteros− Weinstein numbering, a single most conserved residue among the class A GPCRs is designated x.50, where x is the transmembrane helix number. All other residues on that helix are numbered relative to this conserved position) at the tip of TM3, and cysteine of ECL2 in the middle of extracellular loop 2 (ECL2) was created by using SYBYL 6.9.29 The resultant models were subjected to 1000 cycles of steepest descent and 1000 cycles of conjugate gradient energy minimization. The PROCHECK30 program which analyzes the residue-by-residue geometry and overall structure quality was employed to evaluate the stereochemical quality of the energy-minimized models of D1R and D2R. The starting structure of D3R in the simulation study comes from the crystal structure (PDB code: 3PBL). Parameterization of l-SPD and Molecular Docking. The geometry of l-SPD was built based on its crystal structure31 using the SYBYL 6.9 program and optimized at the DFT/ B3LYP/6-31G** level. It features the same configuration and conformation as our previous study.20 The partial atomic charges of l-SPD was determined with the DFT/B3LYP/6- 31G** basis set by using the Gaussian 09 program.32 Parameterization and topology generation of l-SPD for CHARMM force field were performed by using the Paratool plug-in of VMD.33 The automated molecular docking program AutoDock4.0 with AutoDock Tools was employed to probe the possible complex of l-SPD bound with energy-minimized dopamine receptors (D1R, D2R, and D3R). Binding energy prediction was carried out using X-Score.34 All available experimental data in combination with the lowest binding free energy were used to select the initial structure for subsequent MD simulations. System Setup. The orientation of D1R, D2R, and D3R were adjusted using VMD, and their principal axis of symmetry was aligned to the z-axis using Orient plug-in of VMD. Helmut Grubmuller’s SOLVATE program was then used to solvate the proteins followed by elimination of water molecules located in the hydrophobic protein−membrane interface. A pre-equilibrated rectangular patch of POPC membrane bilayer with solvated lipid headgroups was generated by using the Membrane Builder plug-in of VMD. The partially solvated proteins were immersed in the POPC membrane bilayer by matching the center of mass of proteins with that of membrane, respectively. The resulting placement of each receptor in The Journal of Physical Chemistry B Article 8122 dx.doi.org/10.1021/jp3049235 | J. Phys. Chem. B 2012, 116, 8121−8130
The Journal of Physical Chemistry B Article 77A iew(right)of 1-SPD-D3R complex embedded in a hydrated POPC lipid bilayer. nd stick representation(middle). Water molecules are shown as lines, and lipid molecules are re )and phosphorus(orange sphere)atoms in van der Waals representation. Na and CI ions are also All nonpolar hydrogen atoms are omitted for clarity b8, SPDD|R—SPDD2R Time (ns) gure 2. Root-mean-square deviation(rmsd)of the backbone atoms(CA, N, c)of all six simulation systems as a function of simulation time: (a) backbone rmsd for unliganded DRs;( b)backbone rmsd for I-SPD-DRs. membrane was in agreement with the suggestions from the water, ions, and lipid headgroups held fixed to allow the melting Orientations of Proteins in Membranes(OPM) database. 5 of lipid tails. Subsequently, 1000 cycles of minimization Lipids whose phosphorus atom lies within the region of the followed by a 0.5 ns equilibration was effected with harmoni protein and those within 0.6 A of the protein were removed. constraint on protein so as to guide the system to the nearest The protein-membrane complex was further solvated using local energy minimum. During this stage, forces were applied to the VMD Solvate plug-in. Water molecules added inside the water molecules entering the highly hydrophobic membrane ipid bilayer and the proteins were eliminated. The VMD protein interface to prevent this region from hydration. Aft Autoionize plug-in was then employed to create an ionic minimization and equilibration with the protein constrained, concentration of 150 mM NaCl by transmuting random water the harmonic constraint was released to permit the whole brief summary of the configuration of each final system is listed finally conducted for 50 ns for all systeme duction runs were molecules into Nat and Cl and to neutralize the system. a system to equilibrate for 0.5 ns further. Pr in Table Sl, and a snapshot of one of the constructe simulation system(I-SPD bound D3R system) is shown in RESULTS AND DISCUSSION Molecular Dynamics Simulation Parameters. All MD 3D Structures of D1R, D2R, and D3R. Sequence simulations were carried out using the parallel molecular gnment indicated that the sequence identity in trans- dynamics program NAMD 2.6 and the CHARMM27 force membrane helices was 42. 4% between DiR and B2AR and field7with CMAP terms8for all protein molecules, POPC 72.9% between D2R and D3R, respectively(Figure S1).The lipid molecules, and ions along with the TIP3P model for 3D models of dir and D2R were constructed using the newly water molecules. A cutoff of 12 A(switching function starting at resolved agonist-state crystal structure of PAR and antagonist 10 A)for van der Waals interactions was imposed. The particle- state X-ray structure of DR as templates, respectively. The mesh Ewald techni ique was employed to calculate long-range PRoCHECK statistics calculations showed that 100% of the electrostatic forces without cutoff. Periodic boundary con- residues in our DIR and D2R models were either in the most ditions were imposed on all simulations with an integration favored or in the additionally allowed regions of the time step of 2 fs to allow a multiple time stepping algorithms to Ramachandran plot (Figure S2), suggesting that the overall be employed. The simulations were equilibrated as an NPT main-ch d side-chain conformations are much more semble, using the Langevin Nose-Hoover method* to reliable as compared with our previous models. The previous maintain the pressure at 1 atm, and the temperature was kept at models of DiR and d2r were built by taking bovine rhodopsin 310 K by using Langevin dynamics with a very weak friction as the template, which was the only crystallized structure in coefficient GPCR family at that time and had ca. 25% sequence identity in Molecular Dynamics Simulation Protocol. All TM region with DIR and D2R. The overall architecture of our lation systems were first subjected to 1000 cycles of ste present DiR and D2R models is similar to that of our descent and 1000 cycles of conjugate gradient DIR and D2R minimization followed by a 0.5 ns equilibration with models, three intramolecular hydrophobic interaction clusters 81 dx. dolora/o.021/p30492351 Phys. Chem. B2012116,8121-813
membrane was in agreement with the suggestions from the Orientations of Proteins in Membranes (OPM) database.35 Lipids whose phosphorus atom lies within the region of the protein and those within 0.6 Å of the protein were removed. The protein−membrane complex was further solvated using the VMD Solvate plug-in. Water molecules added inside the lipid bilayer and the proteins were eliminated. The VMD Autoionize plug-in was then employed to create an ionic concentration of 150 mM NaCl by transmuting random water molecules into Na+ and Cl− and to neutralize the system. A brief summary of the configuration of each final system is listed in Table S1, and a snapshot of one of the constructed simulation system (l-SPD bound D3R system) is shown in Figure 1. Molecular Dynamics Simulation Parameters. All MD simulations were carried out using the parallel molecular dynamics program NAMD 2.636 and the CHARMM27 force field37 with CMAP terms38 for all protein molecules, POPC lipid molecules, and ions along with the TIP3P model39 for water molecules. A cutoff of 12 Å (switching function starting at 10 Å) for van der Waals interactions was imposed. The particlemesh Ewald technique40 was employed to calculate long-range electrostatic forces without cutoff. Periodic boundary conditions were imposed on all simulations with an integration time step of 2 fs to allow a multiple time stepping algorithms to be employed. The simulations were equilibrated as an NPT ensemble, using the Langevin Nose−́ Hoover method41 to maintain the pressure at 1 atm, and the temperature was kept at 310 K by using Langevin dynamics with a very weak friction coefficient. Molecular Dynamics Simulation Protocol. All simulation systems were first subjected to 1000 cycles of steepest descent and 1000 cycles of conjugate gradient energy minimization followed by a 0.5 ns equilibration with protein, water, ions, and lipid headgroups held fixed to allow the melting of lipid tails. Subsequently, 1000 cycles of minimization followed by a 0.5 ns equilibration was effected with harmonic constraint on protein so as to guide the system to the nearest local energy minimum. During this stage, forces were applied to water molecules entering the highly hydrophobic membrane− protein interface to prevent this region from hydration. After minimization and equilibration with the protein constrained, the harmonic constraint was released to permit the whole system to equilibrate for 0.5 ns further. Production runs were finally conducted for 50 ns for all systems. ■ RESULTS AND DISCUSSION 3D Structures of D1R, D2R, and D3R. Sequence alignment indicated that the sequence identity in transmembrane helices was 42.4% between D1R and β2AR and 72.9% between D2R and D3R, respectively (Figure S1). The 3D models of D1R and D2R were constructed using the newly resolved agonist-state crystal structure of β2AR and antagoniststate X-ray structure of D3R as templates, respectively. The PROCHECK statistics calculations showed that 100% of the residues in our D1R and D2R models were either in the most favored or in the additionally allowed regions of the Ramachandran plot (Figure S2), suggesting that the overall main-chain and side-chain conformations are much more reliable as compared with our previous models. The previous models of D1R and D2R were built by taking bovine rhodopsin as the template, which was the only crystallized structure in GPCR family at that time and had ca. 25% sequence identity in TM region with D1R and D2R. The overall architecture of our present D1R and D2R models is similar to that of our previously reported DRs. Similar to previous D1R and D2R models, three intramolecular hydrophobic interaction clusters Figure 1. Side view (left) and top view (right) of l-SPD-D3R complex embedded in a hydrated POPC lipid bilayer. Protein is shown in cartoon representation, with l-SPD in ball and stick representation (middle). Water molecules are shown as lines, and lipid molecules are represented by sticks with their nitrogen (blue sphere) and phosphorus (orange sphere) atoms in van der Waals representation. Na+ and Cl− ions are also shown as gray and green spheres, respectively. All nonpolar hydrogen atoms are omitted for clarity. Figure 2. Root-mean-square deviation (rmsd) of the backbone atoms (CA, N, C) of all six simulation systems as a function of simulation time: (a) backbone rmsd for unliganded DRs; (b) backbone rmsd for l-SPD-DRs. The Journal of Physical Chemistry B Article 8123 dx.doi.org/10.1021/jp3049235 | J. Phys. Chem. B 2012, 116, 8121−8130
The Journal of Physical Chemistry B Article /-SPD-DIR d DIR /-SPD-DI E2 86420 Time(ns) Time(ns) b 802086420 e D2R SPD- D2R D2R 3 E2一 8。z Time(ns) Time(ns) /-SPD-D3R SPD-D3R e2E 0864 f38号9z Figure 3. Distance fluctuations in TM3-TM6 and n-o during the MD simulations. for both DIR and D2R were found: cluster 1 is surrounded by DiR shortened gradually from ca. 20 A to ca. 12 A at the first TM3, TM6, and TM7; cluster 2 is in TM2 and TM4; and 30 ns and remained stable for the rest 20 ns in its unliganded cluster 3 is in TM3 and TM5. However, differences are state simulation. It indicated that a transition from active state observed at the cytoplasmic face of DiR, with an outward to inactive state is through a continuum of intermediate movement of TMS and TM6 and an inward movement of TM3 conformations; this kind of behavior is in consistency with the and TM7 in the present model. Such movement of trans- postulate proposed by Brian Kobilka et al. that GPCRs are membrane segments is a typical feature of activated B,AR."As more like molecular rheostats that are able to sample the present DIR model was constructed based on the active- continuum of conformations. Upon binding with 1-SPD, the state B2AR, it is presumably in an active state and is more eparation in TM3-TM6 was kept at ca 15 A(above 14 A) uitable for investigating the agonist effect of 1-SPD during the entire simulation(Figure 3a), implying that 1-SPD Cytoplasmic Helical Movement in Agonistic DIR and has the capacity of stabilizing DiR in its active state. In contrast ntagonistic D2R/D3R. The three molecula to the conformational changes induced by 1-SPD in DIR, 1- assembled by embedding DIR, D2R, and D3R in a rectangular SPD-induced conformational changes at the cytoplasmic faces ox comprising a POPC lipid bilayer surrounded by explicit of both D2R and D3R were not remarkable. The separation in water were simulated for 50 ns each. The MD simulations TM3-TM6 in unliganded D2R is relatively stable and is kept at showed that, for all systems, temperature, mass density, and ca. 12 a during the entire simulation. While in I-SPD bound volume were relatively stable after 2 ns. After that, the D2R, the separation in TM3-TM6 fluctuated slightly at around fluctuation scale became much smaller for the rms deviations of 11 a during the first 8 ns, and then it jumped to ca. 12 A. As for he backbone atoms of all simulation systems( Figure 2), D3R, the overall landscape of distance in TM3-TM6 is indicating that the molecular systems behaved well thereafter. stabilized at 10-ll A in both liganded and unliganded D3Rs The most prominent differences in conformation between ( Figure 3b, c). active and inactive GPCRs are resided at the cytoplasmic face, lonic Lock Formation/Breakage in Agonistic D1R and with the outward movement of TM3 and TM6, as exemplified Antagonistic D2R/D3R. It is also established that the basal by the active crystal structures of B2AR-T4L-Nb80 and B2 AR- conformation of GPCRs are maintained by interhelix loops and lexes" and by the active cry noncovalent interactions between side chains of a set of crucial receptor.A systematic examination of currently available residues in TM regions. One of the best characterized GPCR crystal structures revealed that the movement of TM6 examples of noncovalent interaction is a salt bridge formed apart from TM3 is the key event in the process of GPCR between R3. 50, part of the highly conserved(D/E)RY motif in activation(Table S2), as such conformational change opens the transmembrane helix 3, and the conserved E6.30 in trans- docking pocket for G protein. The distances between TM3 and membrane helix 6. This salt bridge network was nicknamed TM6 in the cytoplasmic ends( distance between centroid of the "ionic lock". Though it was originally proposed to account for Ca atoms of four residues located at the tip of each TM)were the activation of B,AR, accumulated studies have demonstrated onitored for both liganded and unliganded DiR, D2R, and that this interaction was crucial in stabilizing the inactive state D3R as shown in Figure 3. The separation in TM3-TM6 for of the rhodopsin family GPCRs. In order to get a more 8124 dx. dolora/o.021/p30492351 Phys. Chem. B2012116,8121-813
for both D1R and D2R were found: cluster 1 is surrounded by TM3, TM6, and TM7; cluster 2 is in TM2 and TM4; and cluster 3 is in TM3 and TM5.20 However, differences are observed at the cytoplasmic face of D1R, with an outward movement of TM5 and TM6 and an inward movement of TM3 and TM7 in the present model. Such movement of transmembrane segments is a typical feature of activated β2AR.4 As the present D1R model was constructed based on the activestate β2AR, it is presumably in an active state and is more suitable for investigating the agonist effect of l-SPD. Cytoplasmic Helical Movement in Agonistic D1R and Antagonistic D2R/D3R. The three molecular systems assembled by embedding D1R, D2R, and D3R in a rectangular box comprising a POPC lipid bilayer surrounded by explicit water were simulated for 50 ns each. The MD simulations showed that, for all systems, temperature, mass density, and volume were relatively stable after 2 ns. After that, the fluctuation scale became much smaller for the rms deviations of the backbone atoms of all simulation systems (Figure 2), indicating that the molecular systems behaved well thereafter. The most prominent differences in conformation between active and inactive GPCRs are resided at the cytoplasmic face, with the outward movement of TM3 and TM6, as exemplified by the active crystal structures of β2AR-T4L-Nb80 and β2ARGs complexes4 and by the active crystal structures of opsin receptor.42 A systematic examination of currently available GPCR crystal structures revealed that the movement of TM6 apart from TM3 is the key event in the process of GPCR activation (Table S2), as such conformational change opens the docking pocket for G protein. The distances between TM3 and TM6 in the cytoplasmic ends (distance between centroid of the Cα atoms of four residues located at the tip of each TM) were monitored for both liganded and unliganded D1R, D2R, and D3R as shown in Figure 3. The separation in TM3-TM6 for D1R shortened gradually from ca. 20 Å to ca. 12 Å at the first 30 ns and remained stable for the rest 20 ns in its unliganded state simulation. It indicated that a transition from active state to inactive state is through a continuum of intermediate conformations; this kind of behavior is in consistency with the postulate proposed by Brian Kobilka et al. that GPCRs are more like molecular rheostats that are able to sample a continuum of conformations.25 Upon binding with l-SPD, the separation in TM3−TM6 was kept at ca. 15 Å (above 14 Å) during the entire simulation (Figure 3a), implying that l-SPD has the capacity of stabilizing D1R in its active state. In contrast to the conformational changes induced by l-SPD in D1R, lSPD-induced conformational changes at the cytoplasmic faces of both D2R and D3R were not remarkable. The separation in TM3−TM6 in unliganded D2R is relatively stable and is kept at ca. 12 Å during the entire simulation. While in l-SPD bound D2R, the separation in TM3−TM6 fluctuated slightly at around 11 Å during the first 8 ns, and then it jumped to ca. 12 Å. As for D3R, the overall landscape of distance in TM3−TM6 is stabilized at 10−11 Å in both liganded and unliganded D3Rs (Figure 3b,c). Ionic Lock Formation/Breakage in Agonistic D1R and Antagonistic D2R/D3R. It is also established that the basal conformation of GPCRs are maintained by interhelix loops and noncovalent interactions between side chains of a set of crucial residues in TM regions.24 One of the best characterized examples of noncovalent interaction is a salt bridge formed between R3.50, part of the highly conserved (D/E) RY motif in transmembrane helix 3, and the conserved E6.30 in transmembrane helix 6. This salt bridge network was nicknamed “ionic lock”. 43 Though it was originally proposed to account for the activation of β2AR, accumulated studies have demonstrated that this interaction was crucial in stabilizing the inactive state of the rhodopsin family GPCRs.43,44 In order to get a more Figure 3. Distance fluctuations in TM3−TM6 and N−O during the MD simulations. The Journal of Physical Chemistry B Article 8124 dx.doi.org/10.1021/jp3049235 | J. Phys. Chem. B 2012, 116, 8121−8130
The Journal of Physical Chemistry B Article and after ligand binding, the separation tor behaviors before triple mechanism of A-SPD. Of the 18 interacting trans- detailed understanding of dopamine re etween the centroid of membrane residues with I-SPD, 9 were conserved in all atoms CZ, NHl, NH2 in R3.50 and that of atoms CD, oEl, complexes and all 18 residues in both D2R and D3R were OE2 in E6.30 termed N-o distance"( Figure S3)was strictly conserved as revealed by sequence alignment in binding onitored and was chosen as an indicator of ionic lock pocket(Figure Sg). Their binding pes were different (Figure Sd-f) Specifically, in D2R and D3R, the large binding The ionic lock formed between R3. 50 in the E/DRY motif pocket 1 located between TM3, TMS, TM6, and TM7 was and E6.30 was proposed to have an essential role in maintaining slightly narrower than the corresponding binding pocket 1 in GPCRs in an inactive state. In the simulation of unliganded DIR; such narrowness was presumably resulted from the D2R, the N-o distance between R3. 50 and E6. 30 was stably relative bulkiness of the imidazole side chain H6.55 in D2R and formed a stable ionic lock with E6. 30 during MD simulation. Additionally, there was also a significant difference in ae onfined below 3. 8 A(Figure 3e); it indicated that R3. 50 D3R with respect to the amide side chain of N6.55 in DIR However, in the case of the 1-SPD-bound D2R, the orientation between Y7. 43 of D2R and D3R and w7. 43 of D1 orresponding N-o distance fluctuated between 3 and 7 A, in which y7. 43 pointed toward extracellular side, while w7.4 indicating that the ionic lock slightly fluctuated upon antagonist was directed toward cytoplasmic side. As a result of such binding. Similar to the unliganded D2R system, N-O distance difference, the binding pocket 2 of both D2R and D3R were of unliganded D3R was well conserved at ca. 3. 8 A within the separated from pocket 1, whereas in DiR, binding pocket 2 was first 38 ns of the simulation, and then it decreased to 3 A after merged with pocket 1. Previous efforts to rationalize the 38 ns, which signified a reinforcement of the ionic lock. While structural basis of selectivity have also focused on binding for the I-SPD bound D3 receptor, the ionic lock ruptured at ca. regions that are not conserved among DR subfamily 26 ns, the corresponding N-O distance even upsurged to sa primary attention being given to ECL2. 45+6 In DIR, C186, (Figure 3f). In terms of unliganded D2R and D3R, the D187, and $188 of ECL2 were found to be in contact with I- important feature in unliganded DiR simulation is that no ionic SPD with their counterparts in D2R and D3R being C182 lock is established( Figure 3d), but the distance in the TM3 1183, 1184 and C181, S182, 1183, respectively TM6 at the intracellular ends decreased to ca 12 A. In contrast Functionally relevant interactions between I-SPD and DIR, the corresponding distance kept almost constant(16-20 A) D2R, and D3R were then investigated to understand the the simulation of I-SPD-bound DIR, showing the important agonistic and antagonistic mechanism of I-SPD. Hydrogen of supports bonding, electrostatic (red dashed lines), and hydrophobic the point that the breakage of ionic lock is one of the features in interactions were observed as major contributions to I-SPD the active dir binding. As a common binding Functionally Relevant Interaction between /-SPD and protonated nitrogen atom to form electrostatic attraction with DIR, D2R/D3R. The predicted binding energies of I-SPD for negatively charged D3. 32 in DIR, D2R, and D3R(Figure 5a- DIR, D2R, and D3R are very well correlated with the c), such a common interaction was in well agreement with the experimental values as shown in Table 1 and Figure 4 fact that D3. 32 plays an essential role in ligand binding to aminergic GPCRs. 45 Table 1. Predicted and Experimental Binding Energies of 1- Despite the common binding features, differences in the SPD with DiR, D2R, and D3R binding interactions of I-SPD with all three receptors were dicted binding experimental binding energy remarkable. In the I-SPD-DiR complex, the binding orientation of I-SPD with DiR was reversed as compared to DIR 11.00 the binding orientations of 1-SPD with D2R and D3R(Figure 60 Sh). Given that the structure of I-SpD has one substructure of dopamine resided on both ends, this reversion of orientation was reasonable (Scheme 1). In addition, this adjustment of binding orientation of I-SPD with DiR was optimal for the formation of two hydrogen bonds(d ring hydroxyl group D OH with side chain of N6.55 and A ring hydroxyl group A-OH with D175 of r2=0.95 interactions were critical for the binding affinity of I-SPD with DIR as experimental results showed that simultaneous substitution of A-OH and d-oh by methoxyl groups resulted in monosubstitution of either A-OH or D-Oh by a methoxyl Experimental binding group resulted in a 6-fold loss of binding affinity(Table S3, (kalmo) entries 1 and 4). This binding mode is further supported by 4.Correlation between predicted binding energy and site-directed mutagenesis in B2AR in which N6.55 was mental binding energy demonstrated to be involved in a hydrogen bonding to the p. Oh group of adrenaline. #-49Furthermore, Manivet and co- Although there are deviations between the experimental value workers showed that N6.55 in the 5-HT2B receptor is involved (-10.80,-9.60, and-1066 kcal/mol)>and predicted data, in direct or indirect 5-HT binding the general trend is that l-SPD binds to DiR and D3R more In I-SPD-D2R and 1-SPD-D3R complexes, the overall strong ly than to D2R interaction landscape is similar to each other but with the DIR, D2R, or of -spd reversed relative to that in dir and d3r were initially compared ng at elucidating the previously mentioned. In both I-SPD-D2R and I-SPD-D3R, the 81 dx. dolora/o.021/p30492351 Phys. Chem. B2012116,8121-813
detailed understanding of dopamine receptor behaviors before and after ligand binding, the separation between the centroid of atoms CZ, NH1, NH2 in R3.50 and that of atoms CD, OE1, OE2 in E6.30 termed “N−O distance” (Figure S3) was monitored and was chosen as an indicator of ionic lock formation. The ionic lock formed between R3.50 in the E/DRY motif and E6.30 was proposed to have an essential role in maintaining GPCRs in an inactive state. In the simulation of unliganded D2R, the N−O distance between R3.50 and E6.30 was stably confined below 3.8 Å (Figure 3e); it indicated that R3.50 formed a stable ionic lock with E6.30 during MD simulation. However, in the case of the l-SPD-bound D2R, the corresponding N−O distance fluctuated between 3 and 7 Å, indicating that the ionic lock slightly fluctuated upon antagonist binding. Similar to the unliganded D2R system, N−O distance of unliganded D3R was well conserved at ca. 3.8 Å within the first 38 ns of the simulation, and then it decreased to 3 Å after 38 ns, which signified a reinforcement of the ionic lock. While for the l-SPD bound D3 receptor, the ionic lock ruptured at ca. 26 ns, the corresponding N−O distance even upsurged to 5 Å (Figure 3f). In terms of unliganded D2R and D3R, the important feature in unliganded D1R simulation is that no ionic lock is established (Figure 3d), but the distance in the TM3− TM6 at the intracellular ends decreased to ca. 12 Å. In contrast, the corresponding distance kept almost constant (16−20 Å) in the simulation of l-SPD-bound D1R, showing the important role of l-SPD in stabilizing the active state of D1R. It supports the point that the breakage of ionic lock is one of the features in the active D1R. Functionally Relevant Interaction between l-SPD and D1R, D2R/D3R. The predicted binding energies of l-SPD for D1R, D2R, and D3R are very well correlated with the experimental values as shown in Table 1 and Figure 4. Although there are deviations between the experimental value (−10.80, −9.60, and −10.66 kcal/mol)15 and predicted data, the general trend is that l-SPD binds to D1R and D3R more strongly than to D2R. The differences in binding pocket lining among D1R, D2R, and D3R were initially compared aiming at elucidating the triple mechanism of l-SPD. Of the 18 interacting transmembrane residues with l-SPD, 9 were conserved in all complexes and all 18 residues in both D2R and D3R were strictly conserved as revealed by sequence alignment in binding pocket (Figure 5g). Their binding pocket shapes were different (Figure 5d−f). Specifically, in D2R and D3R, the large binding pocket 1 located between TM3, TM5, TM6, and TM7 was slightly narrower than the corresponding binding pocket 1 in D1R; such narrowness was presumably resulted from the relative bulkiness of the imidazole side chain H6.55 in D2R and D3R with respect to the amide side chain of N6.55 in D1R. Additionally, there was also a significant difference in the orientation between Y7.43 of D2R and D3R and W7.43 of D1 in which Y7.43 pointed toward extracellular side, while W7.43 was directed toward cytoplasmic side. As a result of such difference, the binding pocket 2 of both D2R and D3R were separated from pocket 1, whereas in D1R, binding pocket 2 was merged with pocket 1. Previous efforts to rationalize the structural basis of selectivity have also focused on binding regions that are not conserved among DR subfamily, with primary attention being given to ECL2.11,45,46 In D1R, C186, D187, and S188 of ECL2 were found to be in contact with lSPD with their counterparts in D2R and D3R being C182, I183, I184 and C181, S182, I183, respectively. Functionally relevant interactions between l-SPD and D1R, D2R, and D3R were then investigated to understand the agonistic and antagonistic mechanism of l-SPD. Hydrogen bonding, electrostatic (red dashed lines), and hydrophobic interactions were observed as major contributions to l-SPD binding. As a common binding feature, l-SPD employed its protonated nitrogen atom to form electrostatic attraction with negatively charged D3.32 in D1R, D2R, and D3R (Figure 5a− c), such a common interaction was in well agreement with the fact that D3.32 plays an essential role in ligand binding to aminergic GPCRs.45 Despite the common binding features, differences in the binding interactions of l-SPD with all three receptors were remarkable. In the l-SPD−D1R complex, the binding orientation of l-SPD with D1R was reversed as compared to the binding orientations of l-SPD with D2R and D3R (Figure 5h). Given that the structure of l-SPD has one substructure of dopamine resided on both ends, this reversion of orientation was reasonable (Scheme 1). In addition, this adjustment of binding orientation of l-SPD with D1R was optimal for the formation of two hydrogen bonds (D ring hydroxyl group D− OH with side chain of N6.55 and A ring hydroxyl group A−OH with D175 of ECL2) (Figure 5a). These two hydrogen-bonding interactions were critical for the binding affinity of l-SPD with D1R as experimental results showed that simultaneous substitution of A−OH and D−OH by methoxyl groups resulted in a ca. 60-fold loss of binding affinity,18 while monosubstitution of either A−OH or D−OH by a methoxyl group resulted in a ∼6-fold loss of binding affinity (Table S3, entries 1 and 4).18 This binding mode is further supported by site-directed mutagenesis in β2AR in which N6.55 was demonstrated to be involved in a hydrogen bonding to the β- OH group of adrenaline.47−49 Furthermore, Manivet and coworkers showed that N6.55 in the 5-HT2B receptor is involved in direct or indirect 5-HT binding.50 In l-SPD-D2R and l-SPD-D3R complexes, the overall interaction landscape is similar to each other but with the orientation of l-SPD reversed relative to that in D1R as previously mentioned. In both l-SPD-D2R and l-SPD-D3R, the Table 1. Predicted and Experimental Binding Energies of lSPD with D1R, D2R, and D3R receptor predicted binding energy (kcal/mol) experimental binding energy (kcal/mol) D1R −11.00 −10.80 D2R −9.98 −9.60 D3R −10.72 −10.66 Figure 4. Correlation between predicted binding energy and experimental binding energy. The Journal of Physical Chemistry B Article 8125 dx.doi.org/10.1021/jp3049235 | J. Phys. Chem. B 2012, 116, 8121−8130