Binding of Capsaicin to the TRPV1 Ion Channel - ACS Publications

Oct 26, 2015 - ABSTRACT: Transient receptor potential (TRP) ion channels constitute a notable family of cation channels involved in the ability of an ...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/molecularpharmaceutics

Binding of Capsaicin to the TRPV1 Ion Channel Leonardo Darré† and Carmen Domene*,†,‡ †

Department of Chemistry, King’s College London, Britannia House, 7 Trinity Street, London SE1 1DB, U.K. Chemistry Research Laboratory, University of Oxford, Mansfield Road, Oxford OX1 3TA, U.K.

‡

S Supporting Information *

ABSTRACT: Transient receptor potential (TRP) ion channels constitute a notable family of cation channels involved in the ability of an organisms to detect noxious mechanical, thermal, and chemical stimuli that give rise to the perception of pain, taste, and changes in temperature. One of the most experimentally studied agonist of TRP channels is capsaicin, which is responsible for the burning sensation produced when chili pepper is in contact with organic tissues. Thus, understanding how this molecule interacts and regulates TRP channels is essential to high impact pharmacological applications, particularly those related to pain treatment. The recent publication of a threedimensional structure of the vanilloid receptor 1 (TRPV1) in the absence and presence of capsaicin from single particle electron cryomicroscopy experiments provides the opportunity to explore these questions at the atomic level. In the present work, molecular docking and unbiased and biased molecular dynamics simulations were employed to generate a structural model of the capsaicin−channel complex. In addition, the standard free energy of binding was estimated using alchemical transformations coupled with conformational, translational, and orientational restraints on the ligand. Key binding modes consistent with previous experimental data are identified, and subtle but essential dynamical features of the binding site are characterized. These observations shed some light into how TRPV1 interacts with capsaicin, and may help to refine design parameters for new TRPV1 antagonists, and potentially guide further developments of TRP channel modulators. KEYWORDS: vanilloid receptor, molecular dynamics simulations, free energy, alchemical transformation, metadynamics

■

INTRODUCTION Transient receptor potential (TRP) channels constitute a superfamily of cation channels characterized by their remarkable diversity in ionic selectivity and activation mechanisms, responding to a wide range of chemical and physical stimuli.1 They participate in a plethora of functions in excitable and nonexcitable cells playing a critical role in sensory physiology in animals ranging from worms to humans.1 Such widespread physiological activity translates into a great diversity in TRP-associated channelopathies involved in an extensive range of human disorders.2 Thus, knowledge of how these proteins respond to physical and chemical stimuli is of direct therapeutic relevance to diseases that affect almost every major organ in the human body. Understanding the structural basis of such response has been hindered by the lack of experimental structural data. However, recent single particle electron cryomicroscopy (cryo-EM) experiments have proven successful in the determination of the three-dimensional structure of the vanilloid receptor 1 (TRPV1) under three conditions: (i) without any agonist,3 (ii) in the presence of the agonist capsaicin, and (iii) in the presence of the agonists resiniferatoxin (RTX, a vanilloid from Euphorbia resinifera) and the spider Double Knot toxin (DkTx).4 TRPV1 is the main representative of a subfamily of thermosensitive TRP channels that enable somatosensory cells to detect temperature changes © 2015 American Chemical Society

in the environment, being activated by pernicious high temperatures. Additionally, tissue damage and inflammation products modulate the channel by decreasing its thermal activation threshold (∼43 °C). As a result of this feature, the TRPV1 channel plays an essential role in the molecular mechanisms responsible for injury-related hyperalgesia and pain hypersensitivity (see Julius et al.5 and references therein). As previously suggested for the family of TRP channels,6 the TRPV1 structure revealed a transmembrane topology that resembles that of voltage-gated channels, and exhibits a 4-fold symmetry arrangement of the subunits around a central permeation pore. Each subunit consists of six transmembrane helices (S1−S6) plus a re-entrant loop−helix domain located between S5 and S6, resulting in a cone shaped arrangement where the selectivity filter (SF) is located. Helices S5 and S6 and the SF constitute the central pore of the channel, which is flanked by a voltage-sensor-like (VS-like) domain formed by the S1−S4 helices.3 The structure also provides details of TRPspecific domains: (i) four ankyrin repeats at the N-terminus, reported to be responsible for multiple ligand binding and Received: Revised: Accepted: Published: 4454

August 20, 2015 October 25, 2015 October 26, 2015 October 26, 2015 DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

Figure 1. Capsaicin and TRPV1 structure. (A) Structure of capsaicin. (B) Structural representation of the TRPV1 channel resolved by cryo-EM4 corresponding to PDB entry 3J5R highlighting the transmembrane domain colored by chain: A (blue), B (red), C (gray), D (orange). The black arrow indicates the binding site of capsaicin at the interface between two monomers (A and B) which is magnified on the right. Four identical capsaicin binding sites are present in TRPV1. (C) Two orientations of capsaicin denoted VRdown and VRup described by Jordt and Julius10 and Gavva et al.11 are depicted at the binding site observed in the cryo-EM structure.4

modulation of the channel sensitivity;7 (ii) the P360−V415 linker, conserved among TRPV subtypes, and which connects the ankyrin repeats to the S1 helix; and (iii) an extended and kinked interfacial helix following S6, which corresponds to the signature “TRP domain” found in many TRP channels.1 Perhaps the most studied agonist of TRPV1 channels is (E)N-[(4-hydroxy-3-methoxyphenyl)methyl]-8-methylnon-6-enamide, commonly known as “capsaicin” (Figure 1A). Capsaicin comes from plants of the genus Capsicum, and it is the active component of chili peppers responsible for the burning sensation. After the discovery that a specific and selective receptor for capsaicin was present in nociceptive neurons,8 and the subsequent characterization of such receptor as the TRPV1 ion channel,9 efforts have been made to understand the molecular basis of the interaction between capsaicin and TRPV1. In particular, the organism-specific sensitivity of TRPV1 to capsaicin has been exploited to identify key residues in the binding process. By swapping domains between rat (sensitive) and chicken (insensitive) TRPV1 proteins, the residue Y511 in the S3 transmembrane helix was found essential for the binding of capsaicin.10 Similar chimeric rat/ rabbit (sensitive/insensitive) constructs led to the identification of a single mutation, I550T, also essential for binding.11 Furthermore, introduction of the mutation T550I was sufficient to substantially reduce the affinity of human as well as rat TRPV1 channels for capsaicin. Based on these mutagenesis experiments, two models of binding have been predicted. In the first one, model VRdown, the vanillyl ring (VR) would point downward, interacting with the aromatic ring of Y511.10 In contrast, in the second model, VRup, the capsaicin molecule would adopt an inverted orientation. The vanillyl ring would point upward and would interact with T550 and W549, and the aliphatic tail would point toward Y511.11 The recently reported structure of TRPV1 solved by cryoelectron microscopy4 describes density assigned to capsaicin located in a pocket involving residues from helices S3 (Y511) and S4 (M547 and T550), in agreement with previous reports on vanilloid binding.10,11 Further contacts were identified including E570 from the S4−S5 linker, and L669 from S6 of an adjacent subunit. However, due to the

resolution of the structure (4.2 Å), the orientation of capsaicin in the binding site was not resolved, precluding any description of its mode of binding. More recently, a combination of docking calculations and functional experiments was used to provide structural insight into the mode of binding of capsaicin, supporting a “tail-up, head-down” conformation, similar to VRdown.12 Such a binding pose would require nonspecific hydrophobic contacts of the aliphatic tail with the channel, including residue F543, and hydrogen-bond interactions between the amide moiety and T550. Additionally, a mechanism for ligand induced activation was proposed, requiring the formation of an extra hydrogen bond between the VR and residue E570.12 Modulation and affinity of capsaicin for the TRPV1 channel is still under debate, particularly due to the channel multimodal responses to environmental stimuli such as heat and protons, which have been shown to affect capsaicin binding.13 Initial measurements of capsaicin binding to TRPV1 channels reported an EC50 of 700 nM,9 corresponding to a relatively low response compared to RTX, for which the EC50 was reported to be 40 nM.9 However, more recent measurements indicate an EC50 of ∼50 nM for capsaicin.14 Additionally, a reduction of one pH unit in the environment has been shown to induce a 5-fold decrease of EC50,13b and 2-fold in the presence of Mg2+.15 In the present work, binding models of capsaicin to TRPV1 were generated using docking and over half a microsecond of accumulated molecular dynamics (MD) simulations based on the new cryo-EM structures. The unbiased sampling obtained from the MD trajectories was enhanced by metadynamics simulations aimed at exhaustively exploring the possible capsaicin poses in the binding pocket and building the corresponding free energy landscape. The picture depicted by these calculations was further used as a starting point for free energy calculations, both free energy perturbation and thermodynamic integration, tailored to estimate the standard free energy of binding.16 The atomic detail obtained from this work constitutes a step forward in the understanding of the interaction between capsaicin and this member of the TRP channel family. In addition to capsaicin, other extracts from plants act as TRPV1 agonists, and although their specific binding sites have 4455

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

molecules. The atomic system was solvated by the Solvate plugin of VMD,21 and water molecules within 1.2 Å of protein, and lipid atoms were removed. Ions were added to the system up to a final concentration of 150 mM NaCl. The CHARMM22/ CMAP force field was used for the protein, CHARMM3622 for the lipids, the TIP3P model23 for water, and the standard CHARMM and NBFIX parameters24 were used for ions.25 Parameters for capsaicin were obtained by dividing the molecule in molecular components and by searching for homologue parameters in the CHARMM General force field (CGenFF)26 (see Supporting Information for details). Initial steric clashes were removed by 5000 steps of minimization followed by a sequence of three 500 ps equilibration runs with constraints applied to relax the channel structure softly in the simulation box. The constraints were applied on (i) the lipids head groups, capsaicin molecules, and all protein atoms; (ii) all protein atoms; and finally (iii) the protein backbone atoms only. MD simulations were performed in the NpT ensemble for 100 ns using NAMD2.9.27 Semi-isotropic pressure coupling at 1 atm was accomplished using the Nosé−Hoover Langevin piston,28 while temperature control was achieved by means of the Langevin thermostat.29 Long-range electrostatic interactions were treated using the particle mesh Ewald algorithm,30 and van der Waals forces were smoothly switched off between 10 and 12 Å. The RATTLE algorithm31 was used to constrain all bonds involving hydrogen atoms, in order to use a 2 fs time step. The multi time step algorithm r-RESPA32 was used to integrate the equations of motion. Nonbonded short-range forces were computed every time step, while long-range electrostatic forces were updated every two time steps. Characterization of the structure of capsaicin in the binding site was accomplished by performing principal component (PC) analysis on the position of the agonist in the structurally aligned binding site during the simulations, followed by hierarchical clustering in the PC space33 by means of the Rbased BIO3D package.34 Metadynamics Simulations. Enhanced sampling MD simulations were performed using well-tempered metadynamics to explore the possible conformations of capsaicin in the TRPV1 binding pocket. In this technique, a history dependent biasing potential acting on a set of collective variables is constructed from the deposition of Gaussians with decreasing height as the simulation proceeds. Application of such bias favors the sampling of regions of the collective variable space difficult to access in normal MD simulations. The free energy surface as a function of the selected collective variables is obtained after completion of a simulation by adding the contribution of the deposited Gaussians. For such purpose, two collective variables were defined in this work: (i) the distance between T550 hydroxyl oxygen and the capsaicin amide nitrogen and (ii) the distance between Y511 hydroxyl oxygen and the capsaicin amide oxygen. Potential energy walls were imposed to both collective variables (upper walls at 10 and 8 Å for collective variables i and ii, respectively) to avoid exploring conformations far from the binding pocket region. In addition, to favor the conformational sampling of the ligand in a binding pocket resembling the experimental structure, potential walls were also imposed on the chi1 torsions of Y511 and T550. NAMD2.927 compiled with the PLUMED 2.02 plug-in35 was used using the same setup and parameters employed in the unbiased MD simulations previously described. The starting point was a snapshot taken from the unbiased MD simulations after equilibration was achieved (∼30 ns). Gaussians with initial

not been identified yet, it is likely that they share a common location with capsaicin. A few examples of members of this group of agonists include piperine (from black pepper), zingerone, gingerol and shogaol (from ginger), and various components found in some essential oils (rose, thyme, palmarosa) such as citronellol and gerianol (see Nilius and Appendino17 and references therein). Accordingly, this study contributes useful information to the understanding of TRP regulation, with paramount pharmacological implications.

■

COMPUTATIONAL METHODS TRPV1 Channel Modeling. The starting coordinates of the TRPV1 channel in the absence (APO) and presence (HOLO) of capsaicin were taken from PDB entries 3J5P and 3J5R, respectively. Only the transmembrane domains (residues: 394− 719) were used in the computational model (Figure 1B) as binding of capsaicin occurs in a pocket located between the voltage-sensor-like domain and the pore, above the S4−S5 linker. Molecular Docking Calculations. Although the HOLO form of the channel solved by EM presents some density assigned to capsaicin, the resolution of the structure hampers the determination of the specific mode of binding. To obtain initial models of the structure of capsaicin in the binding pocket, to be further refined by MD simulations, molecular docking calculations using AutoDock4.018 were performed. A Lamarckian genetic algorithm search and an empirical free energy function were used.19 Docking parameters included an initial population of 200 randomly placed individuals, a maximum number of 3 million energy evaluations, a maximum of 27,000 generations, mutation and crossover rates of 0.02 and 0.80, respectively, and an elitism value of 1. AutoTors (as implemented in Autodock Tool Kit20) was used to define the ligand torsional degrees of freedom. The grid maps representing the protein during the docking process were calculated with AutoGrid 4.0, with grid dimensions of 70 × 70 × 70 grid points, a spacing of 0.375 Å, and the center of the box placed at the reported capsaicin EM density. Thirty independent runs were carried out, and the docking solutions were clustered by root-mean-square deviation (RMSD), using a threshold of 2 Å. The resulting clusters were ranked according to the values of the free energy of binding. The best docked conformations were those found to have the lowest binding energy and the largest number of members in the cluster, indicating good convergence. The docking calculations were performed in one out of the four possible binding sites (one per interface between adjacent monomers). Subsequently, the best docking pose was replicated in all four sites, and the resulting system (VRdownlike system) was refined by MD simulations. For clarity, the following notation for capsaicin molecules and binding sites has been used: sites sE, sF, sG, and sH are occupied respectively by capsaicin molecules CAPE, CAPF, CAPG, and CAPH. These sites correspond to the pocket formed by chains A and B, B and D, C and A, and D and C of the channel respectively (see Figure 1B for a definition). Molecular Dynamics Simulations. Each form of the channel was embedded in a pre-equilibrated lipid bilayer of 1palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) molecules, with the axis aligned to the bilayer normal, and the hydrophobic transmembrane domains aligned with the bilayer. Lipid molecules closer than 1.0 Å to protein atoms were removed, with a resulting bilayer consisting of 490 POPC 4456

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

Figure 2. Free Energy Perturbation approach for the estimation of the free energy of binding of capsaicin to the TRPV1 channel. (A) Thermodynamic cycle used for the calculation. The asterisk indicates that restraints are applied to the ligand L. The gray font denotes decoupling of the ligand from the environment. (B) Definition of the restraints applied to aid in the so-called “wandering ligand” problem and limitations in classical MD sampling. Backbone atoms of residues Y511, M547, and T550 are represented by blue spheres, and the center of mass of these atoms is represented by yellow spheres. Similarly, in the ligand, yellow and red spheres represent the center of mass of selected pairs of atoms and the atoms themselves.

tool38 as implemented in VMD21 was used to analyze the results. To account for microreversibility, the conformation of capsaicin needed to be restrained during the free-energy calculations. For that purpose, a specific conformation was selected from the MD simulations of the HOLO system (see Results and Discussion). In each case, calculations were performed on a capsaicin molecule bound to one of the four sites while the other three remained also occupied. To avoid the “wandering ligand” problem in the bound state, positional and orientational restraints were used as in the work by Gumbart et al.16 Figure 2B shows the definition of the geometrical restraints. Briefly, three reference sites in the protein (P1, P2, P3) were defined corresponding to the center of mass of the backbone atoms of selected residues. In addition, three reference sites in the ligand (L1, L2, L3) were defined corresponding to the center of mass of selected pairs of atoms. These six sites were used to define the harmonic restraints ur (distance P1−L1), uθ (angle P2−P1−L1), and uφ (dihedral P3−P2−P1−L1) to control the position of the ligand in the pocket, and uΘ (angle P1−L1−L2), uΦ (dihedral P1−L1−L2− L3), and uΨ (dihedral P2−P1−L1−L2) to control its orientation. When a harmonic restraint on the ligand RMSD (uc) is also applied, complete control on the capsaicin pose in the pocket is achieved. To estimate the contribution to the free energy of binding as a result of the restraints, thermodynamic integration (TI) simulations were performed coupling the force constant of each harmonic potential to the order parameter λ in a 12-point grid. The gradient of the potential energy with respect to the collective variable was calculated from MD simulations at each value of λ, consisting of 200,000 datacollection steps (0.4 ns), after 50,000 steps (0.1 ns) of equilibration. Scaling of the force constants was also performed bidirectionally, and the free energy contribution of the restraints was retrieved averaging both contributions. The error was estimated as the maximum hysteresis between forward and backward calculations. The complete calculation of the free energy of binding also requires compensating for the positional and orientational restraints applied in the bound form but in the unbound state. For such purpose, direct numerical integration of equations S1 and S2 in the Supporting

height of 0.1 kcal/mol and sigma value of 0.3 Å were deposited every 1 ps, with a bias factor of 10. A total metadynamics simulation time of 340 ns was required. The free energy landscape as a function of the two collective variables biased during the simulation was constructed using the sum_hills utility of PLUMED 2.02. The time-independent free energy estimator proposed by Tiwary and Parrinello36 was used to obtain an estimation of the sampling error associated with the free energy surface and, thus, of the calculation convergence. The error was calculated as the standard deviation between the free energy obtained adding the Gaussians deposited during the well-tempered metadynamics simulation (340,000 Gaussians) and the time-independent estimator. The driver utility of PLUMED 1.3 was used to extract snapshots corresponding to each free energy basin, and the GROMACS g_cluster tool was used to cluster the conformations of capsaicin in each basin. Calculation of Standard Free Energy of Binding. Accurate prediction of the standard free energy of binding of capsaicin to TRPV1 was performed following a protocol described by Gumbart et al.16 This method relies upon alchemical transformations, combined with geometrical restraints required to cope with the problem of incomplete sampling inherent to classical MD simulations. The thermodynamic cycle used in the calculation is shown in Figure 2A. Reversible coupling of capsaicin to its environment, either bulk aqueous solution (unbound state) or the TRPV1 pocket (bound state), was performed using the Free Energy Perturbation (FEP) method, in which the order parameter λ was evenly divided into 45 windows of width equal to 0.02 in the range 0 to 0.9, plus 20 windows of width equal to 0.005, in the range 0.9 to 1.0. Each window consisted of 200,000 datacollection steps (0.4 ns), preceded by 50,000 equilibration steps (0.1 ns) for a total simulated time of 26 ns per coupling/ decoupling cycle. For both the bound and unbound states, the calculation was carried out bidirectionally, creating and annihilating the ligand in separate free-energy calculations. The convergence of the alchemical steps was evaluated by analyzing the overlap of the ΔU distributions for the forward and backward calculations, per lambda window (Figure S3). The statistical data was combined by means of the Bennett acceptance ratio (BAR),37 which provides a maximum-likelihood estimator of the free-energy change. The ParseFEP 4457

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

Figure 3. Comparison between docking, experimental, and MD conformations of capsaicin in the binding pocket of TRPV1. (A) Superposition of the main binding conformers from docking calculations (licorice representation), and the spatial occupation maps (red wire-frame isosurface, occupancy = 0.6; blue wire-frame isosurface, occupancy = 0.1) of capsaicin in sites sE, sF, sG, and sH during the last 70 ns MD simulations. (B) Comparison between the experimental electron density and the occupancy density map from the MD simulation (molecule CAPE).

the geometry of the binding site in the APO form shows a RMSD close to 1 Å, indicating that the pocket samples conformations close to the experimental one without the ligand. An exception is observed for site sE where the magnitude of the RMSD rises to values close to 1.5 Å. This is a consequence of a slight shift of the S4S5 linker toward helix S3, which is related to a chi1 torsion conformational change in Y511 of helix S3 (see discussion below and Figure S7A) that drags I573 present in the S4S5 linker. In the HOLO form, RMSD values of ∼1 Å are observed. Although a slight increase in the RMSD is observed for site sF after ∼70 ns, this change is within the fluctuations observed for the other sites, and it is thus regarded as a consequence of intrinsic thermal fluctuations. These results indicate that the local structure of the binding pocket sampled during the simulations of the APO and HOLO forms resembles their experimental counterparts. The only difference reported between the capsaicin binding site in the APO and HOLO cryo-EM structures was a rotation of the side chain of Y511, proposed to be part of an induced fit binding mechanism.4 To investigate this hypothesis, the orientation of this residue was analyzed with and without capsaicin bound (Figure S7A,B). In general, the orientation observed during the simulation of the APO form corresponds to the one observed in the experimental structure where the lateral chain points away from the binding site (Figure S7C). However, Y511 in site sE alternates between the APO and HOLO experimental orientations during the simulation. This result indicates that spontaneous flipping of the Y511 side chain is possible in the absence of the ligand. In the HOLO system, the side chain of Y511 points toward the binding pocket in all four sites during the entire simulation, in agreement with the experimental structure when capsaicin is bound (Figure S7B,C). These observations suggest a conformational-selection model of binding more than an induced-fit mechanism. Y511 would most probably be in the conformation pointing out of the pocket in the APO form, although it could spontaneously and reversibly flip to the conformation facing the pocket. However, upon ligand binding, the latter conformation would be preferred, stabilizing the ligand in the binding site. To compare the best docking solution with the average conformation of capsaicin during the MD simulation starting from such a binding pose, occupancy maps of the agonist were

Information was performed, and the values obtained were combined with the rest of the contributions in equation S3.16

■

RESULTS AND DISCUSSION Docking Results. Figure S4 shows the distribution of docking results clustered by RMSD and classified by binding energy. Consistent with the VRdown model, the cluster with the lowest energy and highest population presents a vanillyl ring pointing toward Y511 and the S4S5 linker, and an aliphatic tail pointing upward in the direction of F543 in helix S4. This orientation of capsaicin in the binding pocket resembles the “tail-up, head-down” pose recently reported.12 The stability of this docking solution was further confirmed by MD simulations as shown in the following sections. Worth noting is that the second most populated cluster corresponds to the ligand in an inverted orientation compared to the main docking solution, so that the aliphatic tail of capsaicin is pointing to Y511 instead of the vanillyl moiety. This conformation resembles model VRup although the interaction with W549, which was proposed by Gavva et al.,11 is not observed. W549 instead is directed toward the outside of the binding site in the experimental structure. Therefore, although a binding pose consistent with model VRup is conceivable, some of the specific contacts reported by Gavva et al.11 based on a homology model of the channel are not observed in the docking solution based on the experimental structure. In what follows, attention will be focused on the best docking solution, namely, the VRdown model. MD Simulation Results. To evaluate the structural stability of the channel in both the HOLO and APO forms, the RMSD of the Cα atoms from the experimental structures was measured during the simulation (Figure S5). Average values of ∼2.5 Å and ∼3 Å are observed for the APO and HOLO forms, respectively, which are within experimental error (3.2 and 4.2 Å resolution for the APO and HOLO experimental structures, respectively). This is indicative of structural sampling of the transmembrane domain close to the experimental conformation. The stabilization of the structures was achieved after the first 30 ns, which will be considered as the starting point for subsequent analysis. The RMSD of the Cα atoms of the residues constituting each of the four binding sites is shown in Figure S6. In general, 4458

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

Figure 4. Conformational dynamics of capsaicin in binding site sG. PC1 versus PC2 colored by (A) simulation time step (blue, time = 0 ns; red, time = 100 ns) and (B) cluster (black, cluster G1; red, cluster G2). (C) Average capsaicin conformation for cluster G2. Contacts with high permanence times are shown in licorice representation. (D) Residence time matrix for each cluster in site sG. The color scale indicates residence times between 40% and 100%. A contact is considered to occur when the distance between protein and capsaicin heavy atoms is lower than 5 Å. Distribution of ligand conformations mapped on (E) variable V1 (distance Y511OH−capsaicinamideO) and variable V2 (distance T550O−capsaicinamideN) and (F) variable V2 and variable V3 (T550 chi1 torsion).

calculated from the MD trajectory and superimposed to the docking conformer and are shown in Figure 3. For three out of four capsaicin molecules (CAPE, CAPG, and CAPH), the average orientation of the agonist was in agreement with the conformation predicted by docking calculations, as confirmed by the degree of superposition of the latter with the occupancy map (isosurface corresponding to 60% of the simulation time). The similarity between the isosurface shape and the molecular structure of the agonist provides an indication of the fluctuations of capsaicin in the binding site. While for CAPG the isosurface closely resembles the general shape of the molecule, a somewhat less defined shape is observed for CAPE and CAPH, indicating a higher degree of mobility in the binding site. In the case of CAPF, the agonist shows an occupancy map that only roughly resembles the docking pose when an occupancy isosurface corresponding to 10% of the simulation time is used (blue mesh). In contrast, the isosurface of comparable occupancy to the one of CAPE, CAPG, and CAPH (red, isosurface = 60%) reveals a substantial change in

the conformation, with the vanillyl ring moving away from Y511 and the aliphatic tail moving toward the S5 helix. Considering the docking results presented previously, which suggest that the binding pocket might be able to accommodate both VRup and VRdown conformations, the rearrangement observed in CAPF could be indicative of an intermediate state in the switching pathway between these two binding modes. Alternatively, it could indicate that simultaneous occupation of all four binding sites is somewhat disfavored. These partial occupation maps obtained from the MD trajectories were also compared with the electron cryomicroscopy density in Figure 3B. A satisfactory overlap is apparent between the area occupied by capsaicin during the MD simulations and that in the cryo-EM structure. In summary, three sites show an average orientation of capsaicin consistent with the VRdown conformation, while one site shows a substantial rearrangement where the vanillyl ring and the aliphatic tail move in opposite directions to adopt an alternative pose. 4459

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

When comparing CAPE, CAPG, and CAPH, the same orientation of capsaicin is observed, within the fluctuations mentioned above, supporting the stability of the VRdown conformation, as expected from the average behavior revealed by the occupancy maps (Figure 3). Additionally, several of the observed contacts between capsaicin and the channel in sites sE, sG, and sH are in agreement with previous reports including Y511,4,10 M547,4 T550,4,11,12 E570,12 L669,4 and F543.12 One of such reports,12 although it supports the interaction between the amide moiety in capsaicin with T550, it proposes a hydrogen bond with the amide oxygen of capsaicin instead of the nitrogen, as indicated by simulations reported in this work. As previously observed for site sG, long-lived interactions with water molecules inside the pocket are also observed in sites sE and sH, confirming the relevance of such watermediated interactions and building a picture where three main hydration points exist. The first of such sites (hydration site I) is occupied by a water molecule forming simultaneous hydrogen bonds with the amine nitrogen of capsaicin, and residues T550 and N551 of the protein (site sE). The second site (hydration site II) is occupied by a water molecule that bridges the capsaicin amine nitrogen solely with residue T550, and it is localized in the interface between the voltage-sensorlike domain and the pore (sites sH and sG). The third site (hydration site III) is occupied by water molecules stabilized by residues E570 and/or S512; they form hydrogen bonds with the oxygen atom of the methoxy and alcohol groups of the vanillyl ring. The variation of the location of the bridging water molecule interacting with the amine nitrogen of the agonist adds an extra dynamical contact point for capsaicin that also contributes to its local mobility in the pocket. In fact, the presence of these water bridges explains the position of the capsaicin amide nitrogen in the pocket. To quantify such position, the distribution of capsaicin conformations was mapped in site sG (for cluster G2) using two variables: V1 and V2. V1 is the distance between the hydroxyl oxygen in T550 and the amide nitrogen in capsaicin, and V2 is the distance between the hydroxyl oxygen in Y511 and the amide oxygen in capsaicin. Such distribution is depicted in Figure 4E and shows a peak at ∼2.8 Å in V1 and ∼4.5 Å in V2. This is consistent with a hydrogen bond between capsaicin amide oxygen and Y511, and a water-mediated interaction between capsaicin amide nitrogen and T550. A second peak is observed at the same position in V1 but at ∼3.0 Å in V2, indicative of simultaneous hydrogen bonds between capsaicin amide oxygen/nitrogen with Y511/T550, respectively. The same observations are displayed in the case of CAPE and CAPH (Figure S10A,E). However, the peak at 3.0 Å in V2 is missing, indicating lack of a direct hydrogen bond with T550. In addition to the changes in the orientation of Y511 when comparing the APO and HOLO structures during the simulation, a flip in the T550 chi1 torsion is also revealed for the HOLO form. These changes depend on the location of the bridging water molecules and the specific site. The distribution of capsaicin conformations in site sG (for cluster G2) is displayed in Figure 4F as a function of two variables, V2 and V3. V3 is the chi1 torsion of T550. Residue T550 in site sG explores two orientations. In the first one, the T550 methyl group points outside the pocket and the hydroxyl group points to the inside (chi1 ∼ 55) in agreement with the experimental conformation. In the second one, the methyl group points toward the site and the hydroxyl group faces toward the outside of the pocket (chi1 ∼ −55). The first orientation shows two

To gain a deeper understanding about the local rearrangements of the ligand in the binding pocket, the agonist poses were classified by decomposing its spatial fluctuations using principal component (PC) analysis coupled with hierarchical clustering on the PC space. For clarity, only the results for site sG are shown in this article, while the corresponding data for the remaining sites (sE, sF, and sH) is presented in the Supporting Information. The results of the PC-based clustering analysis for capsaicin in binding site sG are summarized in Figure 4. The distribution of ligand conformations during the simulation was mapped on the first two PCs of capsaicin, colored by time step (Figure 4A) or by cluster identity (Figure 4B). Two clusters are observed: one poorly populated (black dots in Figure 4B) and occurring only at the beginning of the simulation (blue dots in Figure 4A), and a second one highly populated (red dots in Figure 4B) and homogeneously sampled (mixture of blue and red dots in Figure 4A). This indicates that capsaicin in site sG samples one main conformational basin with an average structure close to the one depicted in Figure 4C. The supremacy of one conformation is consistent with the recognizable molecular shape observed in the spatial occupation maps shown in Figure 3A. The protein residues close to the ligand were identified by measuring the permanence time of the contacts, as shown in Figure 4D. The main interactions are described in Figure 4C. These include Y511, which appears to form a hydrogen bond with the amide oxygen of capsaicin, and also contributes to the hydrophobic cluster formed by L515, I573, F587, and L669 that anchors the vanillyl ring. A second hydrophobic cluster formed by M547, A665, F591, and F543 interacts with the hydrophobic tail of capsaicin. In addition, long-lived interactions with water molecules in the binding site were also observed (Figure 4C,D). This water molecule forms hydrogen bonds with the amide nitrogen atom of capsaicin, and with the hydroxyl oxygen of T550. Therefore, it seems that specific hydration sites in the binding pocket exist, acting as water bridges between the ligand and residue T550. The PC-based cluster analysis of CAPE and CAPH showed two possible poses each (Figure S8A−C,G−I). In CAPE, the difference between conformations includes a slight shift of the vanillyl ring accompanied by a rotation of the hydroxyl group. Such rotation is accompanied by the loss of an intramolecular hydrogen bond with the capsaicin methoxy group and the formation of an intermolecular hydrogen bond with E570 (Figure S8B,C: transition from the red to the black cluster). A slight rotation of the aliphatic chain is also observed while the tail end remains in the same position. This is in agreement with the more extended occupancy surface observed for CAPE in Figure 3. In the case of CAPH, the change is also subtle, involving a rotation of the vanillyl ring with the methoxy group moving toward Y511, with a rotation in the chi2 torsion of the latter. This is also in agreement with the occupancy map observed for CAPH in Figure 3A, where a narrower vanillyl mesh is indicative of a rotation of this moiety. As indicated by the temporal behavior of the clusters (Figure S8A−C,G−I), capsaicin alternates between clusters during the simulation. Consequently these motions correspond to the intrinsic fluctuations of the agonist in the binding pocket. In fact, although a few specific contacts are formed/lost while sampling these two clusters in CAPE (i.e., F522 or Y554) or CAPH (i.e., S512, N551, or L574), the overall contact pattern is conserved (Figure S9A,C). 4460

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

Figure 5. Free energy landscape of capsaicin in the TRPV1 binding pocket. (A) 2D free energy landscape obtained from well-tempered metadynamics calculations, describing the binding modes of capsaicin in terms of two collective variables: (i) the distance between T550 hydroxyl oxygen and the capsaicin amide nitrogen and (ii) the distance between Y511 hydroxyl oxygen and the capsaicin amide oxygen. (B) Structure of capsaicin in the binding pocket (only helices S3 and S4 and S4S5 linker are shown for clarity) at the four energy minima and the metastable region. Residues Y511 and T550 are depicted in licorice representation. The position of water molecules mediating the interaction of capsaicin with T550 and Y511 is indicated by the black arrows. Hydrogen bonds are indicated with black dotted lines.

interaction peaks with T550−CAPG distances of ∼3.0 Å and ∼4.5 Å, corresponding to direct or water-mediated hydrogen bonds, respectively. In the second orientation, no direct hydrogen bond between CAPG and T550 is observed. In site sE (Figure S10B), the T550 lateral chain remains close to the experimental orientation, with the T550 hydroxyl group pointing to the inside, which is consistent with the occupation of hydration site I. In contrast, in site sH (Figure S10F), the T550 chi1 torsion rapidly shifts to the conformation with the methyl group pointing toward the site and interacting with the

vanillyl ring, while the hydroxyl group faces toward the outside of the pocket. Despite such rotation, a water-mediated hydrogen bond with the amine nitrogen of capsaicin is still possible but with a water molecule localized in hydration site II (Figure S8I). These observations suggest that the orientation of T550 can have an effect on the binding of capsaicin by changing the position of the bridging water molecules, adding a hydrophobic contact with the T550 methyl group and allowing capsaicin to form a direct hydrogen bond with the protein. 4461

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

S14D bottom). Specific localization of water in binding interfaces is known to affect the steric complementarity between protein and ligand, and it mediates binding by providing a hydrogen-bond network.39 Our data provides evidence in favor of such a role of water molecules in the binding of capsaicin to TRPV1. Water accessibility may be also functionally relevant considering that some residues localized in capsaicin binding site (i.e., E570) are susceptible to changes in the pH, which has been shown to regulate gating and capsaicin affinity.13b Free Energy Landscape of Capsaicin in the TRPV1 Binding Pocket. The stability of the VRdown conformation was illustrated by the unbiased MD results described previously. A complex atomic picture was revealed as a consequence of pocket side-chain flexibility and water-mediated hydrogen-bond interactions. Thus, to obtain a more exhaustive and quantitative description, the intrinsically limited sampling of MD simulations was enhanced by employing well-tempered metadynamics calculations. As observed previously, direct or watermediated hydrogen bonds between capsaicin and residues T550 and Y511 constitute pivotal anchoring points. Thus, a twodimensional potential of mean force was built using as collective variables (i) the distance between T550 hydroxyl oxygen and the capsaicin amide nitrogen and (ii) the distance between Y511 hydroxyl oxygen and the capsaicin amide oxygen. Figure 5A shows the free energy landscape describing the localization of capsaicin in the binding pocket of TRPV1 as defined by these two collective variables. Four minima were found and the convergence of their free energy values was accounted for by calculating the sampling error using a timeindependent free energy estimator (see Computational Methods). Minimum 1 corresponds to direct hydrogen bonds between the protein and capsaicin. Minimum 2 corresponds to a conformation of capsaicin that hydrogen bonds with residue Y511, and the distance to T550 is ∼4.5 Å. In minimum 3, capsaicin directly forms a hydrogen bond with residue T550 and is at a distance of ∼5 Å from Y511. Finally, minimum 4 corresponds to capsaicin interacting through a hydrogen bond with residue T550, and at a distance to residue Y511 of ∼7.5 Å. In addition, a metastable region (denoted with 3′ in Figure 5A) is observed, where capsaicin is at distances of ∼4.5 and 5 Å from residues T550 and Y511 respectively. Figure 5B illustrates the binding pocket and the structure of capsaicin in all the identified minima, and in the metastable region; only helices S3 and S4 and S4−S5 linker are shown for clarity. All minima share the same orientation of capsaicin in the binding pocket, consistent with the VRdown conformation, but differences arise in the interactions with residues T550 and Y511. Direct hydrogen bond between capsaicin and residue T550 is observed in three of the minima (1, 3, and 4), while a direct hydrogen bond with Y511 is only observed in two of the minima (1 and 2). In addition, the space between capsaicin and T550 is occupied by a water molecule which bridges the ligand and residue T550 in minimum 2. The same is found in minimum 3, where a water molecule mediates the interaction with Y511. In the case of minimum 4, the wider space between capsaicin and Y511 allows accommodating two water molecules capable of forming a hydrogen-bond network linking Y511 and capsaicin. Finally, the metastable region denoted 3′ corresponds to a binding pose where the interaction with both Y511 and T550 residues occurs via water-mediated hydrogen bonds. These results confirm the importance of hydration sites in the binding of capsaicin to TRPV1, as suggested by unbiased MD

The conformations observed in site sF differ from the others. The PC-based cluster analysis shows four clusters (see Figure S8D−F) that represent the evolution of the ligand from an initial conformation (black representation in Figure S8B) comparable to the other three ligands, to a final conformation where the vanillyl ring has moved to a position above residue T550. The conformation of residue T550 in the first three clusters corresponds to that of the methyl group pointing toward the pocket. However, in the fourth conformation, another shift in T550 chi torsion is observed (Figure S10D), probably as a consequence of the motion of the vanillyl ring. It is worth emphasizing that, although capsaicin undergoes this orientational change in site sF, two of the experimentally observed contacts remain, namely, T550 and L669, indicating their pivotal role in maintaining capsaicin in the binding site. The previous observations depicts an atomic picture of the binding of capsaicin to TRPV1 where subtle changes in the orientation of the residues of the binding pocket and the specific localization of water molecules exert an effect on the local orientation of capsaicin and its contacts with the protein residues, in particular T550. The role of the latter in providing a water-bridged or direct hydrogen-bond contact as observed in the simulations is in agreement with mutagenesis data where the T550I mutant reduces the avidity of the agonist for the channel.10 According to this data, depending on the orientation of T550, at least four scenarios are possible. In the first three, the T550 oxygen is pointing toward the pocket, the experimentally observed orientation, and (i) it collaborates with N551 to anchor a water molecule located in hydration site I, that bridges the protein and the ligand; (ii) it anchors a water molecule between the VS-like and pore domains in hydration site II, that also bridges the protein and the ligand; or (iii) it forms a direct hydrogen bond with capsaicin. In the fourth situation the methyl group of T550 points to the pocket and interacts with the vanillyl ring motif, while the oxygen anchors a water molecule located in hydration site II, although it faces to the outside of the binding pocket. Previously, it was argued that the ∼4.5 Å distance between the amide nitrogen of capsaicin in the pocket with respect to T550 in CAPE, CAPG, and CAPG was the result of hydration sites bridging capsaicin to T550. This was supported by the high permanence times of some water molecules in the neighborhood of capsaicin (Figures 4D and S9). To demonstrate the occurrence of water bridges between ligand and protein, the distance between capsaicin amine nitrogen and T550 hydroxyl oxygen and the distance between the latter two with respect to a water molecule occupying a hydration site were tracked. Figure S14A−C shows such analysis for CAPE, CAPG, and CAPH. For CAPE, the distance between the water and both capsaicin and T550 remain at hydrogen-bond separation during the 30−100 ns of simulation, indicative of a very stable hydration point. In the case of CAPG, initially the agonist interacts directly with T550 at hydrogen-bond distance, but after ∼65 ns, a water molecule enters the site and forms a hydrogen bond with both capsaicin and T550, and the distance between capsaicin and T550 rises to 4.5 Å. In the case of CAPH, a slightly more complex situation is displayed whereby two water molecules interchange their positions at ∼55 ns, and T550 remains at ∼4.5 Å. Water molecules access the binding site by at least two pathways: through the interface formed by (i) helix S4 of one monomer and helix S6 of the adjacent monomer (Figure S14D top) or (ii) helices S3 and S4 from the same monomer (Figure 4462

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics simulations. The most stable pose corresponds to one where both direct hydrogen bonds with T550 and Y511 are present (minimum 1). Substitution of the direct hydrogen bond between capsaicin and residue T550 by a water-mediated interaction (minimum 2) leads to an increase in the free energy of ∼2.5 kcal/mol and ∼1.5 kcal/mol if the water-mediated interaction is with residue Y511 (minimum 3). An increase in the distance between capsaicin and residue Y511 (∼7.5 Å), while keeping the direct hydrogen bond with T550 (minimum 4), leads to the least stable pose, ∼3 kcal/mol higher than the global minima. Errors in the free energy landscape shown in Figure S12 are below 0.3 kcal/mol for the identified minima. The error in the free energy difference between minima is of ∼0.4 kcal/mol, and for the metastable region error values are 0.5−0.75 kcal/mol. In summary, the VRdown conformation of capsaicin in the TRPV1 binding pocket is in fact a combination of four energy basins distinguished by the formation of direct or watermediated hydrogen bonds with Y511 and T550, with free energy differences in the order of the formation/breaking of one hydrogen bond. Standard Free Energy of Binding Calculations. To complement the structural and thermodynamic insight already presented, the standard free energy of binding from the bulk aqueous solution to the binding pocket was estimated. Considering the depth of the pocket in the protein, its position embedded in the lipid membrane, and the lack of information about the entry pathway, a combination of FEP and TI was the chosen approach to estimate the free energy (see Computational Methods). Only knowledge of the conformation of the ligand in the binding pocket is required in this technique. The results presented in the previous section support the existence of the VRdown mode of binding, although it was found to be constituted by multiple minima differing in the interaction mode with residues Y511 and T550. Consequently, to obtain an estimation of the standard free energy of binding for the VRdown conformation, it suffices to target the FEP/TI calculations toward only one of those minima, and integrate such result with the free energy landscape obtained via metadynamics. Therefore, the conformation corresponding to direct hydrogen bond between capsaicin and both T550 and Y511 (minimum 1) was selected, since it is the most stable and thus the more populated. The change in free energy associated with each step of the thermodynamic cycle depicted in Figure 2A and the corresponding estimation of the standard free energy of binding are presented in Table 1. A value of −10.6 ± 1.7 kcal/mol was found for the free energy required to take capsaicin from the aqueous solution and to the TRPV1 pocket in the pose corresponding to minimum 1. Comparison to experimental data is not currently possible as the dissociation constant for capsaicin is not experimentally available. Complete understanding of capsaicin binding to TRPV1 would require the knowledge of the entry pathway of the agonist to the binding pocket, either directly from the solvent or through the cell membrane. The latter hypothesis has recently found support by comparison of MD simulations of the VS-like domain (helices S1−S4) in a lipid bilayer and the translocation potential of mean force of capsaicin across the membrane.40 However, the direct description of the binding pathway remains elusive. Other alternative binding modes could be plausible as suggested by the exploratory docking calculations presented here.

Table 1. FEP/TI Results Decomposed in Contributions as Shown in the Free Energy Cycle Described in Figure 2A Contribution ΔGc,site ΔGΘ,site ΔGΦ,site ΔGΨ,site ΔGθ,site ΔGφ,site ΔGr,site ΔGcouple,site ΔGr,bulk + ΔGa,bulk ΔGo,bulk ΔGdecouple,bulk ΔGc,bulk ΔGbind

VRdown-like ΔG (kcal/mol) −4.5 −0.7 −2.1 −0.5 −0.2 −0.4 −1.5 −38.1 7.0 6.9 10.3 13.2 −10.6

± ± ± ± ± ± ± ±

0.2 0.1 0.7 0.0 0.0 0.0 0.1 0.1

± 0.1 ± 1.5 ± 1.7a

a

The total error is calculated as the square root of the sum of individual squared errors.

■

CONCLUSIONS The binding of capsaicin to the TRPV1 channel has been studied by combining molecular docking, unbiased MD simulations, and free energy methods (metadynamics, Free Energy Perturbation, and Thermodynamic Integration) using the first available structure of a transient receptor potential channel. A comprehensive picture of the thermodynamics governing the binding of capsaicin to the TRPV1 channel for the VRdown conformation has emerged. In good agreement with previous experimental reports, capsaicin was found to bind to a pocket of TRPV1 in a conformation where the vanillyl ring points toward the S4S5 linker while the lipid tail points upward in the direction of the S4 helix. Four free energy basins are found for capsaicin in this orientation; in each of them direct or water-mediated hydrogen bonds are established between the amide moiety of capsaicin and residues Y511 and T550 of the channel. These minima are separated by energy differences comparable to the formation or breaking of a single hydrogen bond. In addition, the standard free energy of binding from the bulk solution to one of these basins was estimated. Unfortunately, none of the TRPV1 antagonists developed until now have progressed beyond phase II trials. However, these calculations shed some light onto how TRPV1 interacts with capsaicin, and may help to refine the design parameters for new TRPV1 antagonists, which will certainly pave the way to further developments in TRP pharmacology.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.molpharmaceut.5b00641. Capsaicin force field parametrization description, free energy calculations, molecular docking results, and MD results (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] Tel: +44-(0)2078487541. Notes

The authors declare no competing financial interest. 4463

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics

■

(20) (a) Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J. AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility. J. Comput. Chem. 2009, 30 (16), 2785−91. (b) Sanner, M. F. Python: A programming language for software integration and development. J. Mol. Graphics Modell. 1999, 17 (1), 57−61. (21) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. J. Mol. Graphics 1996, 14 (1), 33−8. (22) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr.; Pastor, R. W. Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types. J. Phys. Chem. B 2010, 114 (23), 7830−43. (23) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926−935. (24) Beglov, D.; Roux, B. Finite Representation of an Infinite Bulk System - Solvent Boundary Potential for Computer-Simulations. J. Chem. Phys. 1994, 100 (12), 9050−9063. (25) Luo, Y.; Roux, B. Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. J. Phys. Chem. Lett. 2010, 1 (1), 183−189. (26) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31 (4), 671− 90. (27) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26 (16), 1781−802. (28) (a) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. Constant-Pressure Molecular-Dynamics Simulation - the Langevin Piston Method. J. Chem. Phys. 1995, 103 (11), 4613−4621. (b) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101 (5), 4177− 4189. (29) Izaguirre, J. A.; Catarello, D. P.; Wozniak, J. M.; Skeel, R. D. Langevin stabilization of molecular dynamics. J. Chem. Phys. 2001, 114 (5), 2090−2098. (30) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103 (19), 8577−8593. (31) Andersen, H. C. Rattle - a Velocity Version of the Shake Algorithm for Molecular-Dynamics Calculations. J. Comput. Phys. 1983, 52 (1), 24−34. (32) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time-Scale Molecular-Dynamics - Reply. J. Chem. Phys. 1993, 99 (3), 2278−2279. (33) Grant, B. J.; McCammon, J. A.; Caves, L. S.; Cross, R. A. Multivariate analysis of conserved sequence-structure relationships in kinesins: coupling of the active site and a tubulin-binding sub-domain. J. Mol. Biol. 2007, 368 (5), 1231−48. (34) Grant, B. J.; Rodrigues, A. P. C.; ElSawy, K. M.; McCammon, J. A.; Caves, L. S. D. Bio3d: an R package for the comparative analysis of protein structures. Bioinformatics 2006, 22 (21), 2695−2696. (35) (a) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for freeenergy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180 (10), 1961−1972. (b) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185 (2), 604−613. (36) Tiwary, P.; Parrinello, M. A Time-Independent Free Energy Estimator for Metadynamics. J. Phys. Chem. B 2015, 119 (3), 736−742. (37) Bennett, C. H. Efficient Estimation of Free-Energy Differences from Monte-Carlo Data. J. Comput. Phys. 1976, 22 (2), 245−268.

ACKNOWLEDGMENTS The authors acknowledge use of the infrastructure in the Hartree Centre, the EPSRC UK National Service for Computational Chemistry Software (NSCCS), and ARCHER, the UK National Supercomputing Service (http://www.archer. ac.uk).

■

REFERENCES

(1) Venkatachalam, K.; Montell, C. TRP channels. Annu. Rev. Biochem. 2007, 76, 387−417. (2) Nilius, B. TRP channels in disease. Biochim. Biophys. Acta, Mol. Basis Dis. 2007, 1772 (8), 805−12. (3) Liao, M.; Cao, E.; Julius, D.; Cheng, Y. Structure of the TRPV1 ion channel determined by electron cryo-microscopy. Nature 2013, 504 (7478), 107−12. (4) Cao, E.; Liao, M.; Cheng, Y.; Julius, D. TRPV1 structures in distinct conformations reveal activation mechanisms. Nature 2013, 504 (7478), 113−8. (5) Julius, D. TRP channels and pain. Annu. Rev. Cell Dev. Biol. 2013, 29, 355−84. (6) Ramsey, I. S.; Delling, M.; Clapham, D. E. An introduction to TRP channels. Annu. Rev. Physiol. 2006, 68, 619−47. (7) Lishko, P. V.; Procko, E.; Jin, X.; Phelps, C. B.; Gaudet, R. The ankyrin repeats of TRPV1 bind multiple ligands and modulate channel sensitivity. Neuron 2007, 54 (6), 905−18. (8) Szallasi, A.; Blumberg, P. M. Specific binding of resiniferatoxin, an ultrapotent capsaicin analog, by dorsal root ganglion membranes. Brain Res. 1990, 524 (1), 106−11. (9) Caterina, M. J.; Schumacher, M. A.; Tominaga, M.; Rosen, T. A.; Levine, J. D.; Julius, D. The capsaicin receptor: a heat-activated ion channel in the pain pathway. Nature 1997, 389 (6653), 816−24. (10) Jordt, S. E.; Julius, D. Molecular basis for species-specific sensitivity to ″hot″ chili peppers. Cell 2002, 108 (3), 421−30. (11) Gavva, N. R.; Klionsky, L.; Qu, Y.; Shi, L.; Tamir, R.; Edenson, S.; Zhang, T. J.; Viswanadhan, V. N.; Toth, A.; Pearce, L. V.; Vanderah, T. W.; Porreca, F.; Blumberg, P. M.; Lile, J.; Sun, Y.; Wild, K.; Louis, J. C.; Treanor, J. J. Molecular determinants of vanilloid sensitivity in TRPV1. J. Biol. Chem. 2004, 279 (19), 20283−95. (12) Yang, F.; Xiao, X.; Cheng, W.; Yang, W.; Yu, P. L.; Song, Z. Z.; Yarov-Yarovoy, V.; Zheng, J. Structural mechanism underlying capsaicin binding and activation of the TRPV1 ion channel. Nat. Chem. Biol. 2015, 11 (7), 518−24. (13) (a) Neelands, T. R.; Jarvis, M. F.; Han, P.; Faltynek, C. R.; Surowy, C. S. Acidification of rat TRPV1 alters the kinetics of capsaicin responses. Mol. Pain 2005, 1, 28. (b) Ryu, S.; Liu, B.; Qin, F. Low pH potentiates both capsaicin binding and channel gating of VR1 receptors. J. Gen. Physiol. 2003, 122 (1), 45−61. (14) Stein, A. T.; Ufret-Vincenty, C. A.; Hua, L.; Santana, L. F.; Gordon, S. E. Phosphoinositide 3-kinase binds to TRPV1 and mediates NGF-stimulated TRPV1 trafficking to the plasma membrane. J. Gen. Physiol. 2006, 128 (5), 509−22. (15) Ahern, G. P.; Brooks, I. M.; Miyares, R. L.; Wang, X. B. Extracellular cations sensitize and gate capsaicin receptor TRPV1 modulating pain signaling. J. Neurosci. 2005, 25 (21), 5109−16. (16) Gumbart, J. C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? J. Chem. Theory Comput. 2013, 9 (1), 794−802. (17) Nilius, B.; Appendino, G. Spices: the savory and beneficial science of pungency. Rev. Physiol., Biochem. Pharmacol. 2013, 164, 1− 76. (18) Morris, G. M.; Goodsell, D. S.; Halliday, R. S.; Huey, R.; Hart, W. E.; Belew, R. K.; Olson, A. J. Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function. J. Comput. Chem. 1998, 19 (14), 1639−1662. (19) Huey, R.; Morris, G. M.; Olson, A. J.; Goodsell, D. S. A semiempirical free energy force field with charge-based desolvation. J. Comput. Chem. 2007, 28 (6), 1145−1152. 4464

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465

Article

Molecular Pharmaceutics (38) Liu, P.; Dehez, F.; Cai, W. S.; Chipot, C. A Toolkit for the Analysis of Free-Energy Perturbation Calculations. J. Chem. Theory Comput. 2012, 8 (8), 2606−2616. (39) Li, Z.; Lazaridis, T. Water at biomolecular binding interfaces. Phys. Chem. Chem. Phys. 2007, 9 (5), 573−581. (40) Hanson, S. M.; Newstead, S.; Swartz, K. J.; Sansom, M. S. Capsaicin Interaction with TRPV1 Channels in a Lipid Bilayer: Molecular Dynamics Simulation. Biophys. J. 2015, 108 (6), 1425−34.

4465

DOI: 10.1021/acs.molpharmaceut.5b00641 Mol. Pharmaceutics 2015, 12, 4454−4465