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