The transient receptor potential (TRP) channels act as key sensors of various chemical and physical stimuli in eukaryotic cells. Despite years of study, the molecular mechanisms of TRP channel activation remain unclear. To elucidate the structural, dynamic, and energetic basis of gating in TRPV1 (a founding member of the TRPV subfamily), we performed coarse-grained modeling and all-atom molecular dynamics (MD) simulation based on the recently solved high resolution structures of the open and closed form of TRPV1. Our coarse-grained normal mode analysis captures two key modes of collective motions involved in the TRPV1 gating transition, featuring a quaternary twist motion of the transmembrane domains (TMDs) relative to the intracellular domains (ICDs). Our transition pathway modeling predicts a sequence of structural movements that propagate from the ICDs to the TMDs via key interface domains (including the membrane proximal domain and the C-terminal domain), leading to sequential opening of the selectivity filter followed by the lower gate in the channel pore (confirmed by modeling conformational changes induced by the activation of ICDs). The above findings of coarse-grained modeling are robust to perturbation by lipids. Finally, our MD simulation of the ICD identifies key residues that contribute differently to the nonpolar energy of the open and closed state, and these residues are predicted to control the temperature sensitivity of TRPV1 gating. These computational predictions offer new insights to the mechanism for heat activation of TRPV1 gating, and will guide our future electrophysiology and mutagenesis studies.
The mammalian transient receptor potential (TRP) channels are a superfamily of polymodal cation-permeable channels divided into six subfamilies (Montell, 2005) based on sequence homology (TRPC, TRPV, TRPM, TRPP, TRPML, and TRPA). TRP channels are versatile cellular sensors (Clapham, 2003; Voets et al., 2005) activated and regulated by various physical and chemical stimuli such as heat (Caterina et al., 1997, 1999), cold (Story et al., 2003; Karashima et al., 2009; Nilius et al., 2012), voltage (Voets et al., 2004), acid (Tominaga et al., 1998; Jordt et al., 2000), force (Howard and Bechstedt, 2004; Sotomayor et al., 2005; Yin and Kuebler, 2010), exogenous ligands (e.g., capsaicin; Caterina et al., 1997), endogenous agents (e.g., ATP; Lishko et al., 2007), and phosphatidylinositol-4,5-bisphosphate (PIP2; Qin, 2007; Rohacs and Nilius, 2007). TRP channels are among the most actively pursued drug targets (Gunthorpe and Szallasi, 2008; Nilius, 2013) in recent years thanks to their extensive involvements in many intracellular signaling pathways and pathophysiology linked to various diseases (Nilius, 2007; Nilius et al., 2005b).
As a founding member of the TRPV subfamily, TRPV1 forms a homotetramer. Each of its four subunits consists of a six-helix (S1–S6) transmembrane domain (TMD) and an intracellular domain (ICD; see Fig. 1 A). The TMD consists of two functional modules: the S1–S4 voltage-sensing module (Voets et al., 2007) on the channel periphery and the S5–S6 pore module enclosing a central pore (Fig. 1, A and B). The N-terminal part of ICD is an ankyrin repeats domain (ARD; Fig. 1 A) comprised of six ankyrin repeats, which are ubiquitous motifs involved in protein–ligand interactions (Gaudet, 2008a). The C-terminal domain (CTD) of ICD (Fig. 1 A) contains a highly conserved TRP box helix (Montell, 2001), which is implicated in coupling to TRPV1 gating (García-Sanz et al., 2007) and interactions with other proteins/ligands (Numazaki et al., 2003; Prescott and Julius, 2003). At the TMD–ICD interface is a membrane proximal domain (MPD; Fig. 1 A), which was found to harbor a heat-sensing module in the TRPV subfamily (Yao et al., 2011). Alternative heating-sensing sites were found in CTD (Vlachová et al., 2003; Brauchi et al., 2006, 2007) and outer pore (Grandl et al., 2010; Yang et al., 2010; Cui et al., 2012; Kim et al., 2013). It remains uncertain which specific residues are involved in heat activation of TRPV1 and how they communicate with the channel gate to trigger its opening (e.g., via allosteric coupling; Latorre et al., 2007).
To elucidate the molecular mechanisms of TRP channels, tremendous efforts in structural biology (Gaudet, 2008b; Hellmich and Gaudet, 2014) have resulted in low resolution (>13-Å) structures of TRP channels visualized by electron microscopy (Mio et al., 2005, 2007; Maruyama et al., 2007; Moiseenkova-Bell et al., 2008; Shigematsu et al., 2010; Cvetkov et al., 2011; Huynh et al., 2014), and high resolution structures for various fragments of TRP channels solved by x-ray crystallography, including the isolated ARD of TRPV channels (Jin et al., 2006; McCleverty et al., 2006; Lishko et al., 2007; Phelps et al., 2008; Landouré et al., 2010; Inada et al., 2012; Shi et al., 2013). In two landmark papers published in 2013, high resolution (3–4-Å) structures were solved by cryo-electron microscopy (Cao et al., 2013; Liao et al., 2013) for a minimal functional construct of rat TRPV1 in three distinct forms: a closed-channel apo form, an open-channel form bound with DkTx (a peptide toxin) and RTX (a small vanilloid agonist), and an intermediate form bound with capsaicin (an agonist). The apo structure is thought to capture the closed (C) state of TRPV1 channel, and the DkTx–RTX-bound structure likely captures the open (O) state of TRPV1, although other possibilities like a desensitized state cannot yet be ruled out. These structures revealed a dual gate with two constrictions in the channel pore (Fig. 1, A and B): a selectivity filter or upper gate (near residue G643) and a lower gate (near residue I679). Both are closed (open) in the closed-channel (open-channel) structure, whereas the selectivity filter is closed and the lower gate is partially open in the capsaicin-bound structure (Cao et al., 2013; Liao et al., 2013; see Fig. 3 C). Despite great functional insights offered by these structures, to fully understand the activation and regulation mechanisms of TRP channels, we need quantitative knowledge of the underlying dynamics and energetics, which cannot be directly deduced from static structures. It is imperative to obtain and comprehend high resolution structural and dynamic information for key functional states (i.e., C and O state) and the gating transition between them in TRP channels.
Molecular dynamics (MD) simulation is the method of choice for exploring protein dynamics and energetics under physiological conditions with atomic resolution (Karplus and McCammon, 2002). MD was used previously to study various ion channels (Lindahl and Sansom, 2008; Nury et al., 2010, 2011; Zhu and Hummer, 2010; Dong and Zhou, 2011; Delemotte et al., 2012; Du et al., 2012; Vargas et al., 2012; Calimet et al., 2013), including truncated domains (Sotomayor et al., 2005; Susankova et al., 2007) and homology models (Brauchi et al., 2007; Susankova et al., 2007; Pedretti et al., 2011) of TRP channels. However, MD simulation is computationally expensive. Until recently, it was impractical to simulate an ion channel beyond tens of nanoseconds (ns) without using a massively parallelized or special-purpose supercomputer (Jensen et al., 2012). Thanks to rapid progress in computational hardware and software (e.g., use of graphics-processing unit to accelerate MD simulation; Stone et al., 2010), one can now routinely run MD simulation of a large biomolecular system (with >105 atoms) at a speed of >1 ns/day on a single computer node. However, it remains difficult for MD to access microsecond–millisecond time scales relevant to many functionally important conformational transitions of protein complexes (e.g., the gating of a TRP channel).
To overcome the timescale limit of MD simulation, coarse-grained modeling has been vigorously pursued using reduced protein representations (e.g., one bead per amino acid residue) and simplified force fields (Tozzini, 2005, 2010). As a popular example of coarse-grained models, the elastic network model (ENM) is constructed by connecting nearby Cα atoms with harmonic springs (Atilgan et al., 2001; Tama and Sanejouand, 2001; Zheng and Doniach, 2003). Despite its simplicity, the normal mode analysis (NMA) of ENM can be used to predict a few low frequency modes of collective motions, which often compare well with conformational changes observed between experimentally solved protein structures in different functional states (e.g., the gating conformational changes in a pentameric ligand-gated ion channel; Krebs et al., 2002; Taly et al., 2005; Sauguet et al., 2014). Numerous studies have established ENM as a useful and efficient means to probe dynamic mechanisms of protein complexes (including membrane proteins; Bahar et al., 2010) with virtually no limit in timescale or system size (Bahar and Rader, 2005; Tama and Brooks, 2006). ENM has formed the basis for our coarse-grained modeling of various conformational transitions between two given protein structures (Zheng et al., 2007a; Tekpinar and Zheng, 2010, 2013; Zheng, 2010, 2011, 2012; Zheng and Auerbach, 2011; Zheng and Tekpinar, 2012), such as the closed and open conformations of a pentameric ligand-gated ion channel (Zheng and Auerbach, 2011).
In this study, to elucidate the molecular mechanism of TRPV1 gating, we have performed a comprehensive coarse-grained modeling and all-atom MD simulation based on the recently solved high resolution structures of open and closed form of TRPV1 (Cao et al., 2013; Liao et al., 2013). Our coarse-grained NMA has captured two key modes of collective motions involved in the gating transition of TRPV1, featuring a quaternary twist motion of the TMDs relative to the ICDs. Our transition pathway modeling has predicted a sequence of structural motions propagating from the ICDs to the TMDs via the MPDs and the CTDs, leading to sequential opening of the selectivity filter followed by the lower gate in the channel pore (confirmed by our modeling of conformational changes induced by the activation of ICDs). The above findings of coarse-grained modeling are shown to be robust to the perturbation of lipids. Finally, our MD simulation of the ICD has identified key residues that contribute differently to the nonpolar energy in the open and closed state, which are predicted to control the temperature sensitivity of TRPV1 gating. This study is, to our knowledge, the first detailed modeling/simulation study of full-length TRPV1 after the publication of the TRPV1 structures (Cao et al., 2013; Liao et al., 2013). The novel computational predictions from this study will guide our future electrophysiology and mutagenesis studies to investigate the heat-activation mechanism of TRPV1 gating.
MATERIALS AND METHODS
ENM and NMA
In an ENM, a protein structure is represented as a network of Cα atoms of amino acid residues. Harmonic springs link all pairs of Cα atoms within a cutoff distance Rc chosen to be 10 Å (following Hafner and Zheng, 2010). The ENM potential energy is:
where N is the number of Cα atoms, θ(x) is the Heaviside function, dij is the distance between the Cα atom i and j, and is the value of dij as given by a reference structure. The spring constant kij is 1 for nonbonded interactions and 10 for bonded interactions (in arbitrary units).
In NMA, the following eigenproblem is solved for the Hessian matrix H, which is obtained by calculating the second derivatives of ENM potential energy (see Zheng and Auerbach, 2011):
where λm and Wm represent the eigenvalue and eigenvector of mode m. After excluding six zero modes corresponding to three rotations and three translations, we kept 3N-6 nonzero modes, which are numbered from 1 to 3N-6 in order of ascending eigenvalue.
To validate ENM-based NMA, we compared each mode (i.e., mode m) with the observed structural changes Xobs between two superimposed experimental structures by calculating an overlap value Im (Zheng, 2012). Im varies between 0 and 1, with higher value meaning greater similarity. gives the fractional contribution of mode m to Xobs. The sum of over the lowest M modes (named cumulative squared overlap) gives the fractional contribution of these modes to Xobs (Zheng, 2012). To assess the local flexibility at individual residue positions as described by the lowest M modes, we define the following cumulative flexibility:
We performed the above ENM-based NMA using the AD-ENM webserver (http://enm.lobos.nih.gov/start.html).
ENM of TRPV1 in the presence of lipids
We used the CHARMM-GUI PACE CG Builder (Qi et al., 2014) to embed the C-state structure of TRPV1 (Protein Data Bank [PDB] accession no. 3J5P) in a coarse-grained bilayer of POPC lipids and keep 102 POPC molecules within 10 Å of the Cα atoms of TRPV1. Each POPC molecule contains 13 coarse-grained beads as defined in the Martini force field (Marrink et al., 2007). Then we built an ENM of 2,368 Cα beads and 1,326 POPC beads with springs connecting any two beads within a cutoff distance of 10 Å.
Coarse-grained transition pathway modeling by interpolated ENM (iENM)
We previously developed an iENM protocol to construct a transition pathway (i.e., a series of intermediate conformations) between two given protein conformations by solving the saddle points of a double-well potential built from these two conformations (Tekpinar and Zheng, 2010). Here we applied this method to the closed and open conformation of TRPV1 using the iENM webserver (http://enm.lobos.nih.gov/start_ienm.html).
The iENM can be used to model an unknown activated conformation from a known inactivated conformation together with a target conformation of the activation domain. The idea is to progressively transform the activation domain toward its target conformation and let the rest of the protein relax to a series of minimal-energy conformations leading to a final model for the unknown activated conformation. We have recently used this method to build the pre-powerstroke conformation of dynein motor domain from its post-powerstroke conformation after ATP binding (Zheng, 2012). Here, we used this method to model an intermediate conformation of TRPV1 upon activation of the ICDs toward channel opening.
Using the iENM-predicted transition pathway, we can determine the motional order of various protein parts by calculating a reaction coordinate for each part S (denoted RCS; for details see Zheng, 2010). RCS varies from 0 to 1 as the transition proceeds from the beginning to the end conformation. For two protein parts S and S′, if RCS < RCS′ along the pathway, we infer that the motion of S′ precedes that of S. Here, we plotted RCS as a function of fprogress (a fractional progress parameter as defined in Zhu and Hummer, 2009) to track the evolution of RCS along the transition pathway.
Coarse-grained transition pathway modeling of TRPV1 in the presence of lipids
We first used iENM to construct a coarse-grained O-state model of TRPV1 in the presence of lipids. Starting from the C-state model of the TRPV1–lipids system (see above), we used iENM to progressively transform TRPV1 toward the O-state conformation (PDB accession no. 3J5Q) and let the surrounding lipids relax to a minimal-energy conformation. Then we used iENM to construct a transition pathway from the C-state model to the O-state model of the TRPV1–lipids system, and calculated reaction coordinates for various functional parts of TRPV1 (see above). We also used iENM to model an intermediate conformation of the TRPV1–lipids system upon activation of the ICDs toward channel opening (see above).
All-atom MD simulation of ICD
Based on two experimental structures of TRPV1 (PDB accession nos. 3J5P and 3J5Q), we prepared two systems (i.e., C and O state) for the truncated ICD of TRPV1 (residues 111–429 and 690–719). The protonation states of residues with ionizable side chains were determined by the PDB2PQR server (Dolinsky et al., 2007). The hydrogen atoms were added with the visual MD (VMD) program (Humphrey et al., 1996). Both systems were immersed into a rectangular box of water molecules extending 10 Å from the protein in each direction by using VMD. To ensure a physiological ionic concentration of 0.15 M and zero net charge, Na+ and Cl− ions were added to the systems by VMD. The systems were refined with two rounds of energy minimization using the steepest descent method: first, a 5,000-step energy minimization with harmonic constraints (force constant = 10 kcal/mol/Å2) applied to all backbone heavy atoms, then a 5,000-step energy minimization with harmonic constraints (force constant = 1 kcal/mol/Å2) applied to all Cα atoms. The systems were then heated to 300 K over 300 ps by MD with the same constraints as the second minimization. Then the systems were equilibrated for 500 ps with MD performed in the NVT ensemble with the same constraints as in heating. Finally, the systems were subject to production MD runs performed in the NPT ensemble. The Nosé–Hoover method (Martyna et al., 1999) was used with temperature T = 300 K and pressure P = 1 atm. The periodic boundary conditions were applied to the systems. A 10-Å switching distance and a 12-Å cutoff distance were used for nonbonded interactions. The particle mesh Ewald method (Deserno and Holm, 1998) was used to calculate long-range electrostatic interactions. The SHAKE algorithm (Hoover, 1985) was used to constrain the bond lengths of hydrogen-containing bonds, which allowed a time step of 2 fs for MD simulation. The energy minimization and MD simulation were performed with the NAMD program (Phillips et al., 2005) using the CHARMM27 force field (MacKerell et al., 1998) and TIP3P water model (Mackerell et al., 2004). For production MD runs in the C/O state, we applied harmonic constraints (force constant = 1 kcal/mol/Å2) to all Cα atoms to restrict conformational sampling near the experimentally solved backbone structures of TRPV1 (allowing side chains to relax freely), and ran five 20-ns-long MD trajectories (total 100 ns in simulation time).
Energetic analysis of ICD
We calculated the nonbonded energy in C/O state for 5 × 500 snapshots of the ICD during the last 10 ns of five MD trajectories. For each snapshot, we separately calculated the van der Waals (vdW) contribution EvdW based on CHARMM27 force field (MacKerell et al., 1998) and the electrostatic contribution Eelec using the Poisson–Boltzmann method (Gilson and Honig, 1988; Im et al., 1998). We used the CHARMM program (Brooks et al., 2009) to partition EvdW to contributions from individual ICD residues (denoted EvdW,n for residue n, following Li and Zheng, 2013, and Zheng et al., 2013). For the Poisson–Boltzmann calculation, a dielectric constant of 4 was used for the protein interior (Gilson and Honig, 1986; Sharp and Honig, 1990a,b; Olson and Reinke, 2000). A dielectric constant of 80 was used for the exterior aqueous environment. A probe radius of 1.4 Å was used to define the molecular surface corresponding to the dielectric boundary. The salt concentration was set to 0.15 M. All the Poisson–Boltzmann calculations were performed using the PBEQ module (Nina et al., 1997; Roux, 1997; Im et al., 1998) of the CHARMM program (Brooks et al., 1983). The atomic Born radii used here were previously calibrated and optimized to reproduce the electrostatic free energy of the 20 amino acids in MD simulations with explicit water molecules (Nina et al., 1997).
We used the ConSurf web server (http://consurf.tau.ac.il/) to assess the degree of residue conservation within the TRPV subfamily. Using the rat TRPV1 sequence as query, we obtained 194 homologous sequences of TRPV proteins from the UNIREF90 database with CS-BLAST (cutoff E-value = 0.0001; maximum/minimum sequence identity = 95/35%). Based on the multiple sequence alignment of these sequences, each residue position is assigned a grade of 1–9 (see supplemental material for details). If the grade is 8 or 9, then the residue is considered to be well conserved.
Online supplemental material
The supplemental material contains three figures, six videos, two PDB coordinate files for MD trajectories in the C and O state, and results of residue conservation analysis from the ConSurf server.
RESULTS and DISCUSSION
ENM-based NMA predicts key collective motions involved in TRPV1 gating
Starting from the C-state structure of TRPV1 (PDB accession no. 3J5P), we constructed a Cα-only ENM, and then performed NMA, resulting in a spectrum of a total 7,098 normal modes (see Materials and methods). Here, we focus on the lowest 20 modes (for expanded results of the lowest 100 modes, see Fig. S1), each describing a specific pattern of collective motion energetically favored by the given structure (see below). For validation of these modes, we assessed how well they capture the experimentally observed conformational change from the C-state structure to the O-state structure of TRPV1 (PDB accession no. 3J5Q). To this end, we calculated the overlap between each mode and the observed closed-to-open conformational change and the cumulative squared overlap of the lowest 20 (100) modes (see Materials and methods). Encouragingly, the lowest 20 (100) modes have captured 72% (82%) of the observed conformational change, with modes 6 and 13 having the highest overlap (Fig. 2 A). Therefore, as found in many protein complexes (Tama and Sanejouand, 2001; Krebs et al., 2002), the ENM-based NMA offers a good description of the functional conformational change involved in TRPV1 gating.
We also compared the lowest 20 modes with the experimentally observed conformational change from the C-state structure to the capsaicin-bound structure of TRPV1 (PDB accession no. 3J5R). In contrast to the closed-to-open conformational change (see above), only 21% of the capsaicin-induced conformational change was captured by the lowest 20 modes, suggesting that high energy modes are involved in driving this conformational change. Collectively, the above findings imply that different gating transition pathways are induced by capsaicin and DkTx/RTX. Future study will be needed to explore the physiological significance of such differences. Here, we will focus on the closed-to-open conformational transition as induced by DkTx/RTX.
Next, we focused on the two modes that contribute most significantly to the closed-to-open conformational change in TRPV1: Mode 6 (with 0.65 overlap and 42% contribution to the observed conformational change) predicts a quaternary twist motion between the TMDs and the ICDs (i.e., a clockwise rotation of TMDs and a counterclockwise rotation of ICDs in the top view; see Video 1). Based on a dynamic domain partition analysis (Zheng et al., 2007b), the collective motion of this mode can be described in terms of rigid-body rotations between the following dynamic domains (Fig. 1 C):
(a) Four dynamic domains in the TMDs (named T1–T4; Fig. 1 C), with each consisting of the S1–S4 helices of one subunit and the S5–S6 helices of an adjacent subunit (with residue G643 of the selectivity filter at the interface of domains T1–T4; Fig. 1 C).
(b) One joint dynamic domain at the TMD–ICD interface (named I1; Fig. 1 C) comprised of residues 400–415 of the MPDs; residues 690–710 of the CTDs; and the bottom of S1, S5, and S6 helices (with residue I679 of the lower gate located inside domain I1; Fig. 1 C).
(c) Another joint dynamic domain at the TMD–ICD interface (named I2; Fig. 1 C) involving the rest of MPDs and the C-terminal residues 752–762.
(d) Four dynamic domains in the ICDs involving primarily four ARDs (named A1–A4; Fig. 1 C).
The above dynamic domain partition roughly agrees with the functional domain assignment (i.e., ARD, MPD, CTD, S1–S4 helices, and S5–S6 pore helices; Fig. 1, A and B), supporting the functional relevance of the collective motions between these domains described by this mode. Strategically located at the hinge between TMDs and ICDs, domains I1 and I2 may be sensitive to chemical/physical perturbations that affect this twist motion and thereby the gating transition of TRPV1. Interestingly, a similar quaternary twist motion was observed between the TMDs and the extracellular domains of a pentameric ligand-gated ion channel, and was proposed to drive its channel opening (Taly et al., 2005, 2006; Cheng et al., 2006; Zheng and Auerbach, 2011).
Mode 13 (with 0.49 overlap and 24% contribution to the observed conformational change) predicts a counterclockwise rotation of the TMDs and a hinge bending motion of the ICDs (see Video 2). According to the dynamic domain partition analysis (Zheng et al., 2007b), this mode can be described in terms of rigid-body rotations between the following dynamic domains (Fig. 1 D):
(a) Four dynamic domains in the TMDs (named T1–T4; Fig. 1 D), with each consisting of the S1–S4 helices of one subunit and the S5–S6 helices of an adjacent subunit (with residue G643 at the interface of domains T1–T4; Fig. 1 D).
(b) One joint dynamic domain at the TMD–ICD interface (named I1; Fig. 1 D) comprised of residues 310–330 of the ARDs; residues 355–370 of the MPDs; and the bottom of S1, S4, S5, and S6 helices (with residue I679 located inside domain I1; Fig. 1 D).
(c) Another joint dynamic domain at the TMD–ICD interface (named I2; Fig. 1 D) involving the rest of MPDs, the CTDs, and residues 280–310 of the ARDs.
(d) Four dynamic domains in the ICDs involving the rest of four ARDs (named A1–A4; Fig. 1 D).
Compared with mode 6, mode 13 describes a similar twist motion of domains T1–T4 and I1 relative to domain I2 (Fig. 1 D), but a different hinge bending motion of domains A1–A4 (Fig. 1 D). Unlike mode 6, domain I1 of mode 13 includes two disjoint parts, and domain I2 of mode 13 is larger (Fig. 1 D). Together, domains I1 and I2 encompass a larger hinge region in mode 13, involving more residues to sense external perturbations that affect the collective motion of this mode. Similar to mode 6, mode 13 shows no opening in the central pore (see Videos 1 and 2).
Besides the above two modes, other low frequency modes describe more collective motions (for detailed visualization of all the lowest 20 modes, check out the NMA results at the AD-ENM webserver: http://enm.lobos.nih.gov/runoutput/TRPV1/run.output). Together, these modes predict high flexibility toward the N terminus of ARDs, at the N terminus of MPDs, and in parts of S1–S4 helices, and low flexibility in the S5–S6 pore helices near residues G643 and I679 (Fig. 2 C). Consistent with our finding, functional data indicated that the ARD is structurally flexible (Liao et al., 2013) and may undergo large conformational change (Salazar et al., 2008), or even unfold in the absence of ligand (Croy et al., 2004).
In sum, our ENM-based NMA predicted two key modes of collective motions involved in the gating transition from the C state toward the O state. However, both modes show no sign of channel pore opening, which is not favored energetically by the C-state conformation of TRPV1.
NMA of TRPV1–lipid system supports the robustness of TRPV1 dynamics to perturbation of lipids
To assess how lipids around TRPV1 affect its collective motions, we repeated NMA for an ENM of TRPV1 in the presence of lipids (see Materials and methods), and compared it with the results of NMA in the absence of lipids (see above). As expected, the introduction of lipid–TMD interactions increases the energy of collective motions in TRPV1, as indicated by greater eigenvalues of the lowest 20 modes (Fig. 2 B). The local flexibility of TRPV1 is significantly reduced in the TMD (particularly in S1–S4 helices) but changes little in the ICD (Fig. 2 C). Then we calculated the overlap between each mode and the observed closed-to-open conformational change. Encouragingly, the lowest 20 (100) modes have captured 75% (83%) of the observed conformational change, with mode 6 having the highest overlap (Fig. 2 A). Next, we compared the lowest 20 modes solved in the presence and absence of lipids by calculating all-to-all pairwise overlaps between them (Fig. 2 D). Although some modes (such as modes 1–5) are essentially invariant to the perturbation of lipids, some (such as modes 10–12) are strongly perturbed by lipids (Fig. 2 D). Among the two key modes involved in the gating of TRPV1, the old mode 6 (solved in the absence of lipids) is highly similar to the new mode 6 (solved in the presence of lipids) with overlap of 0.84 (Fig. 2 D), whereas mode 13 is less invariant but still relatively robust to lipids (with overlap of 0.62; Fig. 2 D). In sum, despite the perturbation of lipids, the low frequency normal modes remain highly effective in capturing the gating conformational change in TRPV1. In particular, the key modes involved in the gating transition (especially mode 6) are robust to the perturbation of lipids. Our finding is consistent with the notion that the structural dynamics of a protein (including membrane protein) depends, to a large extent, on its native structure. However, our finding does not downplay the key role of a protein’s environment (such as lipids and solvent) in tuning its dynamics.
Transition pathway modeling reveals a sequence of structural events during TRPV1 gating
To determine the structural basis of TRPV1 gating, it is critical to probe the entire conformational transition from the C-state conformation to the O-state conformation (Cao et al., 2013; Liao et al., 2013). To this end, we have used the iENM protocol (see Materials and methods) to generate a transition pathway (consisting of a series of intermediate conformations; see Videos 3 and 4) between the two conformations. To assess the robustness of the predicted pathway, we repeated the C-to-O transition pathway modeling in the presence of more complete CTDs, which yielded similar results (see Fig. S3). By performing a reaction-coordinate (RC) analysis of the predicted transition pathway (see Materials and methods), we determined the motional order of the following numerically labeled functional parts of TRPV1 (Fig. 1, A and B):
(a) ARD (residues 110–357, corresponding to RC1): sensing chemical signals from binding of various ligands such as ATP (Lishko et al., 2007) and calmodulin (Rosenbaum et al., 2004; Lishko et al., 2007);
(b) MPD (residues 358–429, corresponding to RC2): sensing heat that activates TRPV1 gating (Yao et al., 2011);
(d) S5–S6 helices of TMD (residues 560–690, corresponding to RC4): forming the channel pore featuring a dual gate (at residues G643 and I679);
(e) S1–S4 helices of TMD (residues 430–559, corresponding to RC5): corresponding to the voltage-sensing module well studied in a related family of voltage-gated ion channels (Catterall, 2012).
By comparing the RCs of the above parts (RC1–RC5) near the middle of the transition pathway, we found RC1 > RC2 > RC3 > RC4, which implies the following motional order: ARD→MPD→CTD→S5–S6 helices (Fig. 3 A). The early motion of ARD is consistent with the finding of highly flexible ARD by NMA (see above), which allows ligand binding at the ARD to affect early events of the gating transition. The motions of MPD and CTD occur at a very similar pace during the first half of the transition pathway, which is consistent with them belonging to the same dynamic domains as found by NMA (Fig. 1, C and D). The above order is consistent with the general model of TRPV1 gating in which the S5–S6 pore helices open in response to physical/chemical perturbations at the ARD, MPD, or CTD. Interestingly, during the first half of the transition pathway, the motion of S1–S4 helices precedes that of S5–S6 helices (i.e., RC5 > RC4; Fig. 3 A), which features a quaternary twist motion (Video 4) also captured by NMA (Video 1). However, the motion of S1–S4 helices lags behind that of S5–S6 helices later during the transition. Therefore, the twisting S1–S4 helices may drive some of the early motion of S5–S6 helices in the gating transition, which may allow voltage change to activate the opening of TRPV1 channel (Voets et al., 2004; Nilius et al., 2005a).
Next we used the RC analysis to compare the motional order of the selectivity filter (represented by residues G643 from four subunits) versus the lower gate (represented by residues I679 from four subunits). Although both remain closed during the first half of the transition pathway, G643 precedes I679 in opening during the second half of the transition pathway (Fig. 3 A). We note that the capsaicin-bound structure of TRPV1 (featuring a closed selectivity filter and a partially open lower gate; Fig. 3 C) does not obey the above order and is therefore not on the closed-to-open transition pathway.
In sum, our transition pathway modeling has predicted a sequence of structural motions propagating from the ARDs to the TMDs via the MPDs and the CTDs, resulting in sequential opening of the selectivity filter followed by the lower gate.
Transition pathway modeling of TRPV1–lipid system supports the robustness of the gating transition pathway to perturbation of lipids
To assess how lipids affect the gating transition of TRPV1, we performed the transition pathway modeling for TRPV1 in the presence of lipids (see Materials and methods), and compared that with the results obtained in the absence of lipids (see above). Despite moderate changes to the RCs, the introduction of lipid–TMD interactions did not change the motional order ARD→MPD→CTD→G643→I679 as found in the absence of lipids (Fig. 3 B). Among all functional parts, the motion of S1–S4 helices is affected the most by lipids (Fig. 3 B): although leading the S5–S6 helices in the absence of lipids (Fig. 3 A), the S1–S4 helices lag behind the S5–S6 helices during the first half of the transition in the presence of lipids (Fig. 3 B). This may be attributed to a larger drag of the S1–S4 helices than the S5–S6 helices by lipids because the former forms more extensive contact with lipids than the latter. Therefore, our transition pathway modeling of TRPV1 gating is largely insensitive to the perturbation of lipids. This finding is consistent with our previous modeling study of the gating transition pathway of a pentameric ligand-gated ion channel, which found good agreement with the experimental Φ values without including lipids in the modeling (Zheng and Auerbach, 2011). However, our finding does not exclude a likely role of lipids in tuning the gating transition of an ion channel like TRPV1.
Modeling of allosteric conformational changes in TMDs induced by activating conformational changes in ICDs
To probe the structural mechanism of allosteric coupling from the ICDs to the TMDs during TRPV1 gating, we attempted to construct an intermediate structural model of TRPV1 after the activation of ICDs. To this end, we started from the C-state structure and used the iENM protocol (see Materials and methods) to progressively transform the ICDs toward the activated conformation given by the O-state structure, then we let the TMDs relax to a series of minimal-energy conformations (Videos 5 and 6). To enable unhindered opening/closing of the channel pore formed by the S5–S6 helices, we turned off the intersubunit elastic interactions between the S5–S6 helices of the TRPV1 tetramer.
The above modeling predicts a quaternary twist motion of the TMDs relative to the ICDs (Video 6), followed by opening of the selectivity filter (at residue G643) while the lower gate (at residue I679) remains closed (Fig. 3 C). This is consistent with our finding that the selectivity filter opens before the lower gate during the gating transition (see above). Therefore, the activated ICDs are allosterically coupled to the more distant selectivity filter instead of the nearer lower gate, highlighting the key role of collective motions in transmitting intramolecular signals over a large distance. Although the activating conformational changes in ICDs can trigger partial opening of the selectivity filter (upper gate), a complete opening of the channel pore requires additional conformational changes in the rest of TRPV1 (i.e., the TMDs).
To further assess the effect of lipids, we repeated the above modeling for the TRPV1–lipid system (see Materials and methods). Unlike the result for TRPV1 alone (see above), we found that neither gate opens after the activation of ICDs in the presence of lipids (Fig. 3 C). Therefore, the lipids seemed to alter the conformational changes in TMDs induced by the ICD activation. This does not necessarily invalidate our above finding (i.e., opening of the upper gate induced by the activation of ICDs). It is possible that our treatment of the TRPV1–lipid system as a homogeneous ENM with uniform force constant may be inaccurate (e.g., the lipids may be less solid-like and more liquid-like than protein, and protein–lipid interactions may be weaker than intra-protein interactions). Indeed, after removing either lipid–lipid or protein–lipid springs while keeping the short-range repulsive interactions between them to prevent collisions (Tekpinar and Zheng, 2010), we were able to recover the finding of upper-gate opening. Therefore, the ENM description of the TRPV1–lipid system may need to be further refined to accurately describe the allosteric conformational changes in TRPV1 induced by the activation of ICDs.
All-atom MD simulation of ICD in C/O state
To further probe how the heat activation of ICD drives the gating transition in TRPV1, we used all-atom MD simulation to analyze the energetics of ICD in the C/O state. Our goal is to identify hotspot residues in ICD that undergo significant change in energy between the C and O state, which may contribute to the temperature sensitivity of TRPV1 gating. We acknowledge that the contributions from residues outside the ICD will not be considered in this simulation, although they may be important to the temperature sensitivity.
To simulate TRPV1 with limited computing resources, we constructed two all-atom systems for a truncated ICD (in the C and O state) submerged in a box of solvent and ions (see Materials and methods), and then conducted extensive MD simulations (i.e., five 20-ns-long trajectories per system; see Materials and methods). To prevent the ICD from being overly flexible (because of the removal of the rest of TRPV1), we applied positional restraints to the Cα atoms (see Materials and methods) during the MD simulation. The MD trajectories were found to stabilize rapidly (within <10 ns; see Fig. S2 A). The simulation of truncated ICD is warranted because: (a) the ICD is known to sense various signals (including heat) for activation and regulation of TRP channels (Gees et al., 2012); and (b) the TRPV1 structures (Cao et al., 2013; Liao et al., 2013) show only a few contacts between the ICDs of adjacent subunits, and previous studies found the ARD forms a structurally independent functional module (Jin et al., 2006; Zheng, 2013). We acknowledge that the bulk of CTD (after residue 719) is unresolved in the cryo-EM structures of TRPV1 (but could be probed using fluorescence resonance energy transfer measurement; De-la-Rosa et al., 2013), which may contribute additional intra- and intersubunit interactions involving the ICDs.
Energetic analysis of ICD in C/O state
The temperature sensitivity of TRPV1 gating requires a large enthalpy difference between the C and O state (denoted ΔH = HO − HC). To probe the energetic basis of heat activation of TRPV1, we assessed contributions to ΔH from polar (electrostatic) energy and nonpolar vdW energy (see Materials and methods). The two contributions were separately evaluated and averaged based on the last 10 ns of five 20-ns-long MD simulations in the C/O state (see Materials and methods). We repeated the calculation for the last 5 ns of MD trajectories to make sure the results are insensitive to truncation of trajectories.
The difference in total vdW energy between the C and O state is 35.3 kcal/mol per ICD (Table 1), which is statistically significant despite large fluctuations in the vdW energy (see Fig. S2 B). This energy difference is sufficient to account for the large enthalpy difference in TRPV1 gating (ΔH of ∼100 kcal/mol for the tetramer or ∼25 kcal/mol per subunit) as determined from the patch-clamp measurement of temperature dependence of single-channel current (Yao et al., 2010). The contribution to vdW energy difference by the MPD is 18.3 kcal/mol (∼52% of the total vdW energy difference; Table 1). In contrast, the difference in electrostatic energy between the C and O state is only −2.7 kcal/mol (Table 1). Therefore, the temperature sensitivity of TRPV1 gating can be explained by a large change in nonpolar energy within the ICD (especially the MPD), which does not require large conformational change in the ICD (root mean squared deviation of ∼1.3 Å). Limited by the truncated ICD simulation and lack of lipids, we are unable to assess the contributions from TMD, ICD–TMD, intersubunit, and protein–lipid interactions. These contributions will be evaluated in the future based on MD simulations of a full-length TRPV1 in the presence of a lipid bilayer.
To identify key residues of ICD that differentially stabilize the C and O state, we used the CHARMM program (Brooks et al., 2009) to partition the vdW energy to individual residues (see Materials and methods), and calculated the difference (denoted as ΔEvdW,n; Fig. 4 A) between the C and O state for each residue (the distribution of ΔEvdW,n is 0.10 ± 0.496 kcal/mol). Using a 2σ cutoff, we identified the following residues of ICD with significant change in vdW energy between the C and O state (Fig. 4 B):
(a) Residues 344, 355, 370, 406, 408, 409, 410, 412, 413, 415, 694, 695, 709, 713, and 714 have higher vdW energy in the O state than the C state. These residues may contribute significantly to the positive enthalpy difference between these states.
(b) Residues 385, 386, 405, 411, 414, 421, 428, and 696 have lower vdW energy in the O state than the C state. These residues may contribute negatively to the enthalpy difference between these states.
The majority of the above residues are in the MPD, supporting the importance of MPD in controlling the temperature sensitivity of TRPV1 gating as found by one of us (Yao et al., 2011). The three top-contributing residues to the vdW energy difference (T406, R409, and M412) are located at a strategic position interfacing between the MPD, the CTD, and the TMD (Figs. 1 A and 4 B). It will be interesting to investigate the functional importance of these residues with mutagenesis and electrophysiology experiments.
Comparison of predicted residues with literature
Many of the above predicted residues in ICD are well conserved in the TRPV subfamily, including G344, R355, T370, R409, M412, L413, L421, I696, and E709 (see Materials and methods). Some of the predicted residues were at (or in proximity to) the targets of past functional and mutational studies of TRPV1. Gain-of-function mutations S343G/R, I352N/T, I689V, and K710A caused weak or strong toxic phenotypes in yeast (Myers et al., 2008). T370 is known to be phosphorylated by protein kinase A and involved in the sensitization of heat-evoked TRPV1 responses (Winter et al., 2013). Mutants I696A, W697A, and R701A showed a severe effect on voltage- and heat-dependent activation of TRPV1 gating (Valente et al., 2008). The K694A/K698A/K710A triple mutant was found to lose its binding affinity to PIP2 (Grycova et al., 2012). In our future mutagenesis and functional studies, we will focus on those untested residues (especially the ones in MPD). In the future, we will also perform MD simulations of the TRPV1 tetramer (in the presence of solvent, ions, and lipids), which will allow us to improve the energy calculations by accounting for contributions from lipids, TMDs, and intersubunit interactions, etc. Before ending, we would like to point out that the cryo-EM structures used for this computational study have several truncated or unresolved regions (such as the outer-pore turret and part of the CTD), which are potentially important to the TRPV1 channel function. We await a future solution of a more complete TRPV1 structure to computationally probe the roles of these missing parts.
Note: while this manuscript was under review, Darré et al. (2015) published an MD simulation study of the permeation and dynamics of the open pore domain of TRPV1 channel, which nicely complements our MD simulation of the ICD of TRPV1.
This study was funded by grants from American Heart Association (14GRNT18980033) and National Science Foundation (0952736) to W. Zheng, and a National Institutes of Health grant (GM104521) to F. Qin. Computational support was provided by the Center for Computational Research at the University at Buffalo and the Biowulf system at the National Institutes of Health.
The authors declare no competing financial interests.
Kenton J. Swartz served as editor.
ankyrin repeats domain
elastic network model
membrane proximal domain
normal mode analysis
transient receptor potential
van der Waals