Binding of Capsaicin to the TRPV1 Ion Channel - Molecular

Subscriber access provided by NEW YORK MED COLL

Article

Binding of capsaicin to the TRPV1 ion channel Leonardo Darré, and Carmen Domene Mol. Pharmaceutics, Just Accepted Manuscript • DOI: 10.1021/acs.molpharmaceut.5b00641 • Publication Date (Web): 26 Oct 2015 Downloaded from http://pubs.acs.org on October 31, 2015

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 free 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 accessible to all readers and 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.

Molecular Pharmaceutics 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 27

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

Molecular Pharmaceutics

Binding of capsaicin to the TRPV1 ion channel Leonardo Darre1 and Carmen Domene1,2 1

Department of Chemistry, King’s College London, Britannia House, 7 Trinity Street,

London SE1 1DB, UK 2

Chemistry Research Laboratory, Mansfield Road, University of Oxford, Oxford OX1 3TA,

UK

Corresponding author: [email protected] Tel: +44 - (0) 2078487541

Running title Capsaicin-TRPV1 interaction Keywords Vanilloid receptor; Molecular Dynamics Simulations; Free Energy; Alchemical Transformation; Metadynamics

TOC

ACS Paragon Plus Environment


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

ABSTRACT Transient receptor potential (TRP) ion channels constitute a notable family of cation channels involved in the ability of 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 three-dimensional structure of the vanilloid receptor 1 (TRPV1) in absence and presence of capsaicin from single particle electron cryo-microscopy experiments, provides the opportunity to explore these questions at the atomic level. In the present work, molecular docking, 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.


Page 2 of 27

Page 3 of 27

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


INTRODUCTION Transient receptor potential (TRP) channels constitute a super family 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 non-excitable 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 threedimensional 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 thermo sensitive TRP channels that enable somatosensory cells to detect temperature changes 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 al5 and references therein). As previously suggested for the family of TRP channels,6 the TRPV1 structure revealed a trans-membrane topology that resembles that of voltage-gated channels, and exhibits a four-fold symmetry arrangement of the subunits around a central permeation pore. Each subunit consists of six trans-membrane 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, S6 and the SF constitute the central pore of the channel, which is flanked by a voltage-sensor-like (VSlike) 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 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,



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

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-3methoxyphenyl)methyl]-8-methylnon-6-enamide, or 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.

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 coloured by chain: A (blue), B (red), C (grey), D (orange) The black arrow indicates the binding site of capsaicin at the interface between two monomers (A and B) which is zoomed in 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


Page 4 of 27

Page 5 of 27

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


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 cryo-electron 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 binding pose would require non-specific 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 has 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 five-fold decrease of EC5013b, and two-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 starting point for free energy calculations, both free energy perturbation and thermodynamic integration, tailored to estimate the standard free energy of binding.16 The



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

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 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 27000 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 70x70x70 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


Page 6 of 27

Page 7 of 27

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


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 resulted system (VRdown-like 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 1-palmitoyl-2oleoyl-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 than 1.0 Å were removed, with a resulting bilayer consisting of 490 POPC molecules. The atomic system was solvated by the Solvate plug-in 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-12 Å. The RATTLE algorithm31 was used to constrain all bonds involving hydrogen atoms, in order to use a 2 fs time-step. The



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

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 R-based 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 set-up 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 height of 0.1 Kcal/mol and sigma value of 0.3 Å where 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 Parrinello,36 was used to obtain an estimation of the sampling error associated to 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


Page 8 of 27

Page 9 of 27

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


during the well-tempered metadynamics simulation (340 000 Gaussians) and the timeindependent 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 to1.0. Each window consisted of 200,000 data-collection 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 bi-directionally, 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 tool38 as implemented in VMD21 was used to analyze the results. To account for micro-reversibility, 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 section). 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 centre 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 centre 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φ



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

(dihedral P3-P2-P1-L1) to control the position of the ligand in the pocket, and uΘ (angle P1L1-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 in 200,000 data-collection steps (0.4 ns), after 50,000 steps (0.1ns) of equilibration. Scaling of the force constants was also performed bi-directionally, 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 Information was performed, and the values obtained were combined with the rest of the contributions in equation e S3.16

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 grey font denotes decouple 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.

RESULTS AND DISCUSSION


Page 10 of 27

Page 11 of 27

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


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 al11 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, 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 raises 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



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

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 resemble their experimental counterparts. The only difference reported between the capsaicin binding site in the APO and HOLO cryoEM 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 alternate 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 binding pose, occupancy maps of the agonist were 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


Page 12 of 27

Page 13 of 27

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


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 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 cryo-microscopy density in Figure 3.B. It is apparent a satisfactory overlap 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.

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).

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



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

principal component (PC) analysis coupled with hierarchical clustering on the PC space. For clarity, only the results for site sG are shown in the main text, 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 intra-molecular hydrogen bond with the capsaicin methoxy group and the formation of an inter-molecular 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.


Page 14 of 27

Page 15 of 27

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


Figure 4. Conformational dynamics of capsaicin in binding site sG. PC1 versus PC2 coloured by (A) simulation time-step (blue: time = 0 ns, red: time=100ns) 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 colour-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 – Capsaicinamide O) and variable V2 (distance T550O-Capsaicinamide N) and (F) variable V2 and variable V3 (T550 Chi1 torsion).

This is also in agreement with the occupancy map observed for CAPH in Figure 3.A where a narrower vanillyl mesh is indicative of a rotation of this moiety. As indicated by the temporal



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

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). 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 L6694 and F543.12 One of such reports,12 although 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 water-mediated 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-sensor-like 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 inY511 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


Page 16 of 27

Page 17 of 27

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


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 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. 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



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

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, as well as 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


Page 18 of 27

Page 19 of 27

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


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 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-chains 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 water mediated hydrogen-bonds between capsaicin and residues T550 and Y511 constitute pivotal anchoring points. Thus, a two dimensional 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 time-independent free energy estimator (see Methods). Minimum 1 corresponds to direct hydrogen-bonds between the protein and capsaicin. Minimum 2 corresponds to a conformation of capsaicin that hydrogenbonds 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 at a distance to residue Y511 of ~ 7.5 Å. In addition, a meta-stable 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.



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

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


Page 20 of 27

Page 21 of 27

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


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, S4 and S4S5-linker are shown for clarity) at the four energy minima and the meta-stable 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.

Figure 5B illustrates the binding pocket and the structure of capsaicin in all the identified minima, and in the meta-stable region; only helices S3, 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 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 meta-stable 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 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 meta-stable 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 water-mediated



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

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 Methods section). 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 to 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. Contribution

VRdown-like ∆G (Kcal/mol)


Page 22 of 27

Page 23 of 27

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


∆Gc,site

-4.5 ± 0.2

∆GΘ,site

-0.7 ± 0.1

∆GΦ,site

-2.1 ± 0.7

∆GΨ,site

-0.5 ± 0.0

∆Gθ,site

-0.2 ± 0.0

∆Gφ,site

-0.4 ± 0.0

∆Gr,site

-1.5 ± 0.1

∆Gcouple,site

-38.1 ± 0.1

∆Gr,bulk+∆Ga,bulk

7.0

∆Go,bulk

6.9

∆Gdecouple,bulk

10.3 ± 0.1

∆Gc,bulk

13.2 ± 1.5

∆Gbind

-10.6 ± 1.7a

Table 1. FEP/TI results decomposed in contributions as shown in the free energy cycle described in Figure 2A. aThe 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 into how TRPV1 interacts with capsaicin, and may help to refine



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

the design parameters for new TRPV1 antagonists which will certainly pave the way to further developments in TRP pharmacology. 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). SUPPORTING INFORMATION 1. Capsaicin force field parameterisation description (Figures S1 and S2. Tables S1 and S2). 2.

Details of Free Energy Calculations (Figure S3).

3. Molecular docking results (Figure S4) 4. MD supplementary results: channels RMSD (Figure S5); binding sites RMSD (Figure S6); Y511 chi1 torsion dynamics in the APO and HOLO forms (Figure S7); PC-based hierarchical cluster analysis for sites sE, sF and sH (Figure S8); Contact permanence time matrices per PC-based cluster in sites sE, sF and sH (Figure S9); Y511/T550capsaicin distance and T550 Chi torsion for sites sE, sF and sH (Figure S10); Hydration sites in the capsaicin-binding pocket (Figure S11); Free energy landscape error estimation (Figure S12).

References 1. Venkatachalam, K.; Montell, C., TRP channels. Annu Rev Biochem 2007, 76, 387417. 2. Nilius, B., TRP channels in disease. Biochimica et biophysica acta 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. Annual review of cell and developmental biology 2013, 29, 355-84. 6. Ramsey, I. S.; Delling, M.; Clapham, D. E., An introduction to TRP channels. Annual review of physiology 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.


Page 24 of 27

Page 25 of 27

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


8. Szallasi, A.; Blumberg, P. M., Specific binding of resiniferatoxin, an ultrapotent capsaicin analog, by dorsal root ganglion membranes. Brain research 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. The Journal of biological chemistry 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. Nature Chemical Biology 2015, 11 (7), 518-+. 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. Molecular pain 2005, 1, 28; (b) Ryu, S.; Liu, B.; Qin, F., Low pH potentiates both capsaicin binding and channel gating of VR1 receptors. The Journal of general physiology 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. The Journal of general physiology 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. Journal of Neuroscience 2005, 25 (21), 5109-5116. 16. Gumbart, J. C.; Roux, B.; Chipot, C., Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? Journal of Chemical Theory and Computation 2013, 9 (1), 794-802. 17. Nilius, B.; Appendino, G., Spices: the savory and beneficial science of pungency. Reviews of physiology, biochemistry and pharmacology 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. Journal of computational chemistry 1998, 19 (14), 16391662. 19. Huey, R.; Morris, G. M.; Olson, A. J.; Goodsell, D. S., A semiempirical free energy force field with charge-based desolvation. Journal of computational chemistry 2007, 28 (6), 1145-1152. 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. Journal of computational chemistry 2009, 30 (16), 2785-91; (b) Sanner, M. F., Python: A programming language for software integration and development. Journal of molecular graphics & modelling 1999, 17 (1), 57-61. 21. Humphrey, W.; Dalke, A.; Schulten, K., VMD: visual molecular dynamics. Journal of molecular graphics 1996, 14 (1), 33-8, 27-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. The journal of physical chemistry. B 2010, 114 (23), 7830-43.



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

23. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of Simple Potential Functions for Simulating Liquid Water. Journal of Chemical Physics 1983, 79 (2), 926-935. 24. Beglov, D.; Roux, B., Finite Representation of an Infinite Bulk System - Solvent Boundary Potential for Computer-Simulations. Journal of Chemical Physics 1994, 100 (12), 9050-9063. 25. Luo, Y.; Roux, B., Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. Journal of Physical Chemistry Letters 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-690. 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. Journal of computational chemistry 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. Journal of Chemical Physics 1995, 103 (11), 4613-4621; (b) Martyna, G. J.; Tobias, D. J.; Klein, M. L., ConstantPressure Molecular-Dynamics Algorithms. Journal of Chemical Physics 1994, 101 (5), 4177-4189. 29. Izaguirre, J. A.; Catarello, D. P.; Wozniak, J. M.; Skeel, R. D., Langevin stabilization of molecular dynamics. Journal of Chemical Physics 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. Journal of Chemical Physics 1995, 103 (19), 85778593. 31. Andersen, H. C., Rattle - a Velocity Version of the Shake Algorithm for MolecularDynamics 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. Journal of Chemical Physics 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 free-energy 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. Journal of Physical Chemistry 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. 38. Liu, P.; Dehez, F.; Cai, W. S.; Chipot, C., A Toolkit for the Analysis of Free-Energy Perturbation Calculations. Journal of Chemical Theory and Computation 2012, 8 (8), 2606-2616.


Page 26 of 27

Page 27 of 27

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


39. Li, Z.; Lazaridis, T., Water at biomolecular binding interfaces. Physical Chemistry Chemical Physics 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. 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.


Binding of Capsaicin to the TRPV1 Ion Channel - Molecular

Recommend Documents