The Molecular Mechanism Underlying Ligand Binding to the

Apr 16, 2018 - ... to the Membrane-Embedded Site of a G-Protein-Coupled Receptor .... Rapid Sampling of Hydrogen Bond Networks for Computational ...
1 downloads 3 Views 990KB Size
Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Biomolecular Systems

The Molecular Mechanism Underlying Ligand Binding to the Membrane-Embedded Site of a G-Protein-Coupled Receptor Xiao-Jing Yuan, Stefano Raniolo, Vittorio Limongelli, and Yechun Xu J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00046 • Publication Date (Web): 16 Apr 2018 Downloaded from http://pubs.acs.org on April 21, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1

The Molecular Mechanism Underlying Ligand

2

Binding to the Membrane-Embedded Site of a G-

3

Protein-Coupled Receptor

4

Xiaojing Yuan†‡, Stefano Raniolo§, Vittorio Limongelli§¶, and Yechun Xu†*

5

†CAS Key Laboratory of Receptor Research, Drug Discovery and Design Center, Shanghai

6

Institute of Materia Medica, Chinese Academy of Sciences (CAS), Shanghai 201203, China

7

‡University of Chinese Academy of Sciences, Beijing 100049, China

8

§Faculty of Biomedical Sciences, Institute of Computational Science - Center for Computational

9

Medicine in Cardiology, Università della Svizzera Italiana (USI), CH-6900 Lugano, Switzerland

10

¶Department of Pharmacy, University of Naples “Federico II”, I-80131 Naples, Italy

11

*Corresponding Author, E-mail: [email protected]

12

13

ABSTRACT The crystal structure of P2Y1 receptor (P2Y1R), a class A GPCR, revealed a

14

special extra-helical site for its antagonist, BPTU, which locates in between the membrane and

15

the protein. However, due to the limitation of crystallization experiments, the membrane was

16

mimicked by use of detergents and the information related to the binding of BPTU to the

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

1

receptor in the membrane environment is rather limited. In the present work, we conducted a

2

total of ~7.5 µs all-atom simulations in explicit solvent using conventional molecular dynamics

3

and multiple enhanced sampling methods, with models of BPTU and a POPC bilayer, both in the

4

absence and presence of P2Y1R. Our simulations revealed that BPTU prefers partitioning into

5

the interface of polar/lipophilic region of the lipid bilayer before associating with the receptor.

6

Then, it interacts with the second extracellular loop of the receptor and reaches the binding site

7

through the lipid-receptor interface. In addition, by use of funnel-metadynamics simulations

8

which efficiently enhance the sampling of bound and unbound states, we provide a statistically

9

accurate description of the underlying binding free energy landscape. The calculated absolute

10

ligand-receptor binding affinity is in excellent agreement with the experimental data (∆Gb0_theo =

11

-11.5 kcal mol−1, ∆Gb0_exp= -11.7 kcal mol−1). Our study broadens the view of the current

12

experimental/theoretical models and our understanding of the protein-ligand recognition

13

mechanism in the lipid environment. The strategy used in this work is potentially applicable to

14

investigate ligands association/dissociation with other membrane-embedded sites, allowing

15

identification of compounds targeting membrane receptors of pharmacological interest.

16 17

INTRODUCTION

18

G protein-coupled receptors (GPCRs) constitute the largest and the most important membrane

19

protein family, and ~40% of marketed drugs act by binding to GPCRs.1 Despite the structural

20

and functional diversity of GPCRs, the ligand binding pockets on the transmembrane region have

21

mostly been found inside the helical bundle.2 However, in 2015, the first extra-helical binding

22

site was reported for an antagonist of the P2Y1 receptor (P2Y1R), BPTU (Figure 1).3 Soon after

ACS Paragon Plus Environment

2

Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1

the discovery of the site in this class A GPCR, such allosteric binding sites outside the helical

2

bundles have also been found in class B GPCRs such as the glucagon-like peptide-1 receptor

3

(GLP-1R) and the glucagon receptor (GCGR).4-6 The possible universal existence of such a

4

binding site among multiple GPCRs makes it intriguing and necessary to study the binding

5

mechanism of ligands with these sites. Since they are located in-between receptors and the

6

membrane, it is likely that the membrane plays a crucial role in protein-ligand recognition

7

processes. However, incompatible with the importance, in crystallization experiments the

8

membrane was usually mimicked by various detergents (such as monoolein) which have no ionic

9

headgroups, unlike the major component of the natural membranes, phospholipids such as POPC

10

lipids (Figure 1). This raises the question if the ligands such as BPTU can still find their ways to

11

bind to the extra-helical sites with the same poses as they were in the experimentally determined

12

structures, when the lipids of the natural environment are used.

13

Moreover, in comparison with the common solvent-exposed sites to which the binding of

14

ligands is governed by the desolvation as well as interactions with receptors,7-10 when the

15

pockets locate inside the membrane, the ligand binding event could be significantly influenced

16

by the lipids and becomes more complicated. Unfortunately, such an important pharmaceutical

17

process is poorly understood in the present time. Several models have been proposed to explain

18

the membrane related phenomena in drug binding. For instance, a “microkinetic” model suggests

19

that the membrane acts as a repository to accumulate drugs around the membrane-embedded

20

receptors, while a “reduction in dimensionality” effect presumes that drugs move from 3D

21

diffusion in the solvent to 2D diffusion within the plane of the membrane so as to increase the

22

drug-receptor collision rate, and the rebinding effect may increase the in vivo duration of the

23

drug effect.11-14 In the context of the microkinetic model, the concentration of hydrophobic

ACS Paragon Plus Environment

3

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

1

ligands in the membrane solvated layer may be driven up due to their affinity to the membrane,

2

with the increase of ligand occupancy in the target compartment. In this way, the experimentally

3

measured affinity may be more potent than the “actual” affinity to the receptor. This membrane

4

accumulation effect thus could muddy the structure-affinity relationship (SAR) analysis in drug

5

design.15 However, these theoretical models are difficult to validate directly with general

6

experimental methods. Although the membrane association effect may be mimicked by use of

7

the immobilized artificial membrane (IAM) method,16 the experiment was conducted without

8

receptors and the direct connection between the membrane effect and the receptor could not be

9

established. Apart from the accumulation effect of the membrane, attentions should also be paid

10

on more advanced questions as to how the membrane could affect drug’s diffusion and what kind

11

of pathway the drug utilizes to associate with its pocket on the receptor.

12 13

Figure 1. Crystal structure of the P2Y1R-BPTU complex, and chemical structures of BPTU,

14

monoolein and POPC. The protein is represented by yellow cartoon and transparent molecular

15

surface, and BPTU is shown in blue sticks. The grey ovals indicate the upper and lower

16

boundary of the membrane.

ACS Paragon Plus Environment

4

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Journal of Chemical Theory and Computation

Due to the difficulty in measuring ligand binding interactions with membrane proteins by use

2

of experimental methods, molecular dynamics (MD) simulations serve as an alternative

3

approach. Previously, long-timescale conventional molecular dynamics (CMD) simulations and

4

enhanced-sampling techniques have been used to study the binding of ligands to the orthosteric

5

pockets of GPCRs,7-8, 17-20 in which the ligands all progressed from the water phase into the

6

extracellular crevice within the helical bundles and the lipids did not play fundamental roles in

7

the binding process. To date, only few researches have investigated the association of ligands,

8

including cholesterol, with GPCRs via the membrane environment.21-26 However, the binding

9

sites in these simulations are mostly intra-helical, the intermediate states along the binding path

10

could be different from those of the extra-helical sites, and the significant role of the lipids has

11

not been clearly characterized yet. A comprehensive interpretation of a ligand binding

12

mechanism needs a detailed description of the interplay among water, lipids, the ligand and its

13

receptor, and an accurate characterization of the underlying free energy landscape. Noteworthy,

14

the access of a ligand to its binding site occurs in the timescale which is usually beyond the

15

capability of CMD, more sophisticated techniques are thereby required. Metadynamics is able to

16

enhance sampling and reconstruct the free energy surface (FES) afterwards, thus it has been

17

widely used in ligand binding studies.8, 18, 27-28 With this method, a history-dependent bias

18

potential that acts on a selected number of degrees of freedom, called collective variables (CVs),

19

is introduced to discourage the system from revisiting the sampled configurations. The well-

20

tempered metadynamics (WT-MetaD) is an evolution of standard metadynamics in which the

21

height of the bias potential introduced is decreased with the amount of bias already deposited.

22

Consequently, the resulting FES could be limited to the low free energy regions which are

23

physically meaningful.29 However, an accurate estimation of the binding affinity needs several

ACS Paragon Plus Environment

5

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32

1

recrossing events between the bound and unbound states, which remains challenging to WT-

2

MetaD simulations due to the large amount of unbound conformations to be sampled. A recently

3

developed metadynamics-based approach, called funnel metadynamics (FM), has proven to be

4

successful in describing ligand/protein and ligand/DNA binding interactions. FM places a

5

funnel-shaped restraining potential on the system so as to enclose the bound states in the cone

6

region and restrain the sampling of the solvated states in a cylinder. As such, a number of back-

7

and-forth events between the bound and unbound states within an affordable computational time

8

become available, which results in a description of a statistically accurate FES and a calculation

9

of the accurate binding affinity.30-32

10

P2Y1R, a class A GPCR, is activated by ADP to induce platelet activation and is thus a

11

promising target for design of novel antithrombotic drugs. BPTU, 1-(2-(2-(tert-butyl)

12

phenoxy)pyridin-3-yl)-3-(4-(trifluoromethoxy)phenyl)urea, is a novel antagonist of P2Y1R and

13

was once developed for the treatment of thrombosis.33 Due to the unique binding site of BPTU

14

and the importance of P2Y1R as a drug target, this system has been recently studied with various

15

computational methods. In 2016, a CMD study was done in order to interpret how BPTU acts as

16

an antagonist of the receptor.34 Meanwhile, based on the BPTU-P2Y1R structure, Xu et al. used

17

the 2D structural similarity search, pharmacophore based screening, and molecular docking to

18

identify new antagonists of P2Y1R from Traditional Chinese Medicine (TCM).35 Later,

19

Antonella et al. performed molecular docking together with CMD studies on P2Y1R and several

20

antagonists, including BPTU, to demystify the P2Y1R-ligand recognition.36 Nevertheless, in

21

these studies BPTU are always positioned in the binding site, the comprehensive binding process

22

of BPTU from a fully solvated state to its final pose in the crystal structure of the complex has

23

not been investigated yet. To clearly reveal the molecular mechanisms underlying the BPTU-

ACS Paragon Plus Environment

6

Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1

P2Y1R recognition, fundamental questions such as how BPTU interacts with water/lipids

2

molecules before approaching the receptor, how BPTU initially associates with the receptor, how

3

BPTU could pave its way to the binding site in the complicated water-lipid environment, how

4

BPTU interacts with key residues to lie itself into the pocket with a favorable binding pose, and

5

what the exact role of lipids are in the recognition process, must be addressed.

6

By combining CMD and multiple enhanced-sampling approaches including umbrella

7

sampling, well-tempered metadynamics and funnel-metadynamics simulations, we explored the

8

molecular mechanism underlying the binding of this allosteric antagonist to the extra-helical site

9

of P2Y1R. Remarkably, We compared the performance of BPTU in the POPC bilayer with or

10

without the receptor, and explicitly discussed the role of lipids in the binding pathway. With

11

these valuable trajectories, we successfully captured the whole association event and

12

characterized the underlying free energy landscape. The current study broadens our

13

understanding of the pharmaceutically relevant ligand/receptor recognition process, especially in

14

the case of binding between ligands and the lipid-exposed sites of membrane proteins, which

15

remains challenging to most experimental techniques.

16

METHODS

17

Simulations of BPTU trafficking inside the POPC bilayer in the absence of the receptor

18

CMD Simulations. The simulation system was built via CHARMM-GUI server by randomly

19

placing one BPTU in aqueous solution with a distance of ~30 Å from the ligand to the center of

20

the bilayer.37 The simulation box, with a size of ~49 Å × 49 Å × 94 Å, contains 72 POPC

21

molecules, 3,875 TIP3P water molecules, and NaCl at a concentration of 0.15 M, resulting in a

22

total of 21,292 atoms approximately. The lipids, water molecules and ions were modeled using

ACS Paragon Plus Environment

7

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

1

the CHARMM36 force field.38 The parameter of BPTU was assigned with the CHARMM

2

General Force Field (CGenFF) by analogy through the ParamChem service (https://cgenff.

3

paramchem.org).39 The system was first minimized and equilibrated using GROMACS40 v5.1.2

4

in a NPT ensemble at 310 K and 1 bar using the Berendsen coupling, and the harmonic-position

5

restraints were applied to the heavy atoms of BPTU and POPC. These restraints were tapered off

6

gradually over the equilibration process. Periodic boundary conditions were applied to the

7

simulation systems. Bonds containing hydrogen atoms were restrained with the LINCS algorithm

8

to allow an integration time step of 2 fs.41 Electrostatic interactions were calculated using the

9

particle-mesh Ewald (PME) summation algorithm.42 In the production runs, the Parrinello-

10

Rahman coupling and the Nose-Hoover coupling were used to maintain a constant pressure of 1

11

bar and a constant temperature of 310 K. Finally, two 500-ns production runs, starting with

12

different initial velocities assigned on each atom of the simulation systems, were resulted.

13

Potential of Mean Force Simulations. The free energy profile of transferring BPTU from the

14

center of the POPC bilayer to the water phase was calculated using the umbrella restraint

15

potential of mean force method. Since pulling from membrane center outward has been shown to

16

increase convergence of final free energy profile,43 one BPTU molecule was placed at the center

17

of the bilayer at the beginning. The size of the box and the following energy minimization and

18

equilibration steps are the same as in the CMD simulations, except that the system was

19

equilibrated for 5 ns in the NPT ensemble. BPTU was then pulled from z = 0 Å (at the center of

20

the bilayer) out to 39 Å (in the bulk water phase) with a pulling rate of 1 Å/ns and a force

21

constant of 1000 kJ/mol/nm2. Snapshots were then extracted from the pulling simulation with

22

BPTU positioned at z = 0 Å, 1 Å, 2 Å, ..., 39 Å from the bilayer center. Following the

23

10ns/window equilibration, 80ns/window production runs were performed with an umbrella

ACS Paragon Plus Environment

8

Page 9 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

1

potential using a force constant of 1000 kJ/mol/nm2 applied to the centroid of the BPTU

2

molecule using the pull geometry cylinder method. The free energy profile was then constructed

3

using the weighted histogram analysis (WHAM) method implemented in Gromacs 5.1.2 with a

4

tolerance of 1 × 10−8.44 The resulted profile was symmetrized and statistical errors were

5

estimated using the bootstrap analysis (Figure S2). As shown in Figure S1, a window of 80ns is

6

long enough for simulations to converge.

7

H-Bonds Analysis. The average numbers of H-bonds between BPTU and water/lipids were

8

calculated for each window of the umbrella sampling using hbond implemented in Gromacs

9

5.1.2, and the default criteria (distance