Transient receptor potential (TRP) ion channels exist in excitable and non-excitable cells and respond to various stimuli, e.g., capsaicin, toxins, heat and pH changes.^[1] Mammal TRP channels can be divided into 6 subfamilies, TRPA, TRPC, TRPM, TRPML, TRPP, and TRPV.^[2] These TRP channels constitute the key sensory elements of human to ambient environment. Mutations of TRP channels can cause lots of human diseases such as autosomal dominant segmental glomerulosclerosis,^[3,4] hypomagnesemia and hypocalcemia.^[5,6] Similar to the K_V channel, TRP channel forms tetramer and has six-transmembrane helixes (S1–S6).^[7] TRP channels also have some unique structural elements such as TRP box which adjuncts to S6 helix and a large ankyrin repeats at the intracellular N terminus.^[8,9] Among the TRP family channels, a group of TRP channels can be activated by heat and form the fundamental sensors of mammals to ambient temperature. These TRP channels include TRPV1-4, TRPM8, TRPM3-5, TRPC5, and TRPA1.^[10–12]

TRPV1 channel (the founding member of TRPV subfamily) is an important member of the thermal-activation channels and mainly exists in dorsal root ganglia neurons and skin cells.^[13,14] TRPV1 has very high Q₁₀ value (∼ 26),^[15] which reflects its high thermal sensitivity. Other chemical stimuli, such as capsaicin,^[13,16] toxin,^[17,18] H⁺, endocannabinoid lipid, and even divalent cation,^[19,20] can also activate TRPV1 channels.^[2] Because TRPV1 channels play an important role in pathological pain,^[21–23] it attracts extensive investigations in the past 20 years. However, the activating mechanism of TRPV1 channels is still controversial, especially the heat-activation mechanism.

TRPV1 channels can be activated by moderate heat (≥ 43 0).^[13] Considering that heat can affect the whole protein and lipid system, it is difficult to identify a specific site for temperature to influence. However, great efforts are taken to verification of the structural basis of the thermal sensitivity from experimental works. From mutagenesis study, the outer and inner pore regions of TRPV1 are confirmed to be important for heat activation.^[24–27] The C-terminal domain is also considered to be important for heat activation because exchanging this part between TRPV1 and TRPM8 (which is a cold-activation channel^[28]) can reverse their hot and cold sensitivity.^[14,28] The turret region which is highly variant among TRP family is mentioned to be part of the temperature-activation pathway.^[8] Chimeric analysis from Yao et al.^[29] revealed that the thermal sensitivity of TRPV1 may arise from a specific protein domain rather than integration of global thermal effects and they confirmed that part of the N-terminal region is crucial for temperature sensing.

Due to the difficulties of direct measurement of the energy for TRPV1’s temperature gating in experiments,^[30,31] different theoretical methods are used to investigate the conformational changes of the secondary structures in heat-activation process and to explain the experimental data. By fitting the experimental data to the two-state model, Voets et al.^[32] showed that the opening rate of TRPV1 is steeply temperature-dependent, whereas the closing rate shows shallow temperature dependence. The kinetic and energetic analyses of Yao et al.^[30] on their direct, time-resolved measurements of TRPV1 showed the energy landscape of gating by heat. The emergence of the high resolution cryo-EM structures of TRPV1 in different states is the milestone of the research of TRPV1 channel.^[33,34] Based on these structures, extensive molecular dynamics (MD) simulations are performed to investigate the activation mechanism by heat or by agonist. Zheng et al.,^[35] using a coarse-grained model, showed a twist motion of the transmembrane domains relative to the intracellular domains in the open process of TRPV1. Their all-atomic simulations of different structures at different temperatures revealed more details of conformational changes in TRPV1’s activation by heat.^[36,37] Simulations with truncated structures of TRPV1 revealed that the large temperature-dependent conformational transition of TRPV1 in the pore domain is asymmetric, only one monomer with large movement can open the pore.^[38] MD simulations are also used to investigate the access of capsaicin to its binding pocket through membrane.^[39] Although the mechanism of TRPV1 gating by heat is much clearer now, a full description of the mechanical pathway from initial conformational change of temperature sensor to the opening of channel gate is still lacking.

We utilized the high resolution cryo-EM structures of TRPV1 to investigate its heat-activation mechanism by using molecular dynamics simulations. By changing the simulation temperatures, we obtained the temperature-dependent opening of the channel pore. The upper gate (G643) suggested in the cryo-EM structures is in the open state in all the simulations (even at 310 K which is under the threshold of heat-activation temperature). M644 has multi-state conformations. Its side chain can occasionally bend to the pore region and block the ion permeation. Opening of the lower gate is temperature-dependent. The temperature sensor of TRPV1 is near the capsaicin-binding site. The temperature-dependent forming and breaking of two hydrogen bonds (HBs) between E570 and S510/S512 induce the oscillation of S5 helix. Through the hydrophobic contact between S5 and S6 helix, the oscillation movement of S5 can transmit to S6 and induces the same movement of S6. Outward movement of S6 finally opens the lower gate. We propose that this is the mechanical pathway from temperature sensor of TRPV1 to the channel opening.

The atomic model of TRPV1 (PDB ID: 3J5P^[34]) which represents a minimal functional construct of rat TRPV1 in the closed state obtained from electron cryo-microscopy experiment is used for modeling. The missing residues (L503-F507) of the original structure are added in the modeling. The highly divergent pore turret (N604-S626) is also added in the calculation model since it is believed to be part of the temperature activation pathway.^[8] The membrane building tool of VMD^[40] is used to generate the bilayer of palmitoyloleoyl-phosphatidylcholine (POPC) and the tetrameric structure of TRPV1 is embedded in this POPC bilayer. The protein-bilayer system is immersed into a box of water molecules with size 98 Å × 98 Å × 96 Å. TIP3P^[41] is used to model water molecules. Na⁺ and Cl⁺ ions are added to ensure an ionic concentration of 150 mM and zero net charge. The entire model has ∼ 350000 atoms. During the simulation run, the α-carbons of residues in S1, S2, S3 and S4 helixes (V430-R455, R474-R499, S512-S532 and V538-Y554) are harmonically constrained because they maintain stationary in channel-activation process.^[33] The software used for modeling and data analysis is VMD (version 1.9.3).^[40] The figures are produced by VMD and Origin.

The MD simulations are performed using NAMD (version 2.12 with CUDA)^[42] with force field CHARMM36^[43–46] at constant temperatures of 310 K, 340 K, 360 K, and 410 K, respectively. The integration time step is 2 fs. The Particle Mesh Ewald (PME) method is used to calculate full electrostatics. Constant pressure control is achieved using Nosé–Hoover Langevin piston at 1 atm. The van der Waals interactions are calculated with a cutoff using a switch function starting at distance of 10 Å and reaching 0 at 12 Å. Langevin dynamics and Langevin piston are performed to maintain the constant temperature and pressure. We performed 10000 steps of energy minimization prior to MD run to remove bad contact. Before the production run, 750 ps simulations with everything except lipid tails is fixed, 750 ps simulations with protein constrained and 750 ps simulations with protein released are performed to equilibrium the system. Simulations in different temperatures are all performed three times with durations of 30 ns, 30 ns and 100 ns. All trajectories are recorded every 2 ps.

To elucidate the open event of TRPV1’s upper and lower gate, the states of the two gates of each simulation frame are calculated. The van der Waals radius of Ca²⁺ is used as the criterion in the calculation because the radius of Ca²⁺ is larger than Cl⁻, and Na⁺, K⁺, and Ca²⁺ can cross the TRPV1 channel. In the CHARMM36^[43–46] force field, the van der Walls radius of Ca²⁺ is 1.367 Å. For the upper gate, the constriction is formed by two O atoms (the radius is 1.7 Å) of G643 at the opposite side. Thus, if the distances between both O–O pairs of opposite G643 are larger than 6.134 Å (1.7 × 2 + 1.367 × 2), we output 1 to denote that this gate is in the open state. For the lower gate, if the distances between the two methyl groups of I679 and that of the opposite I679 are larger than 6.854 Å (2.06 × 2 + 1.367 × 2), the lower gate denotes being in the open state.

The principle component analysis (PCA) of S5 helix is performed using the VMD plugin “Normal Mode Wizard”.^[47,48]

Because E570 is a polar residue with a charged side chain, the criterion for the formation of the hydrogen bond is that the maximum distance between acceptor and donor atoms of the hydrogen bond is 3 Å.^[49] The formation rates of both the hydrogen bonds are indicated by the ratio of the formation time versus the total simulation time (including one 100 ns simulation trajectory and two 30 ns simulation trajectories).

The correlation of outward movement of S5 and S6 helixes is represented by the movement of two C_α atoms of D576 (S5 helix) and I679 (S6 helix). Because we only consider the outward movement of the two C_α atoms which contribute to the lower-gate opening, the displacement vector of the two C_α atoms was projected to the outward movement vector in our simulation. Then the two projected displacements of D576 and I679 are used to calculate the Pearson Correlation Coefficient using the equation

$$ \begin{eqnarray}C=\displaystyle \frac{\displaystyle {\sum }_{i=1}^{n}({X}_{i}-\bar{X})({Y}_{i}-\bar{Y})}{\sqrt{\displaystyle {\sum }_{i=1}^{n}{({X}_{i}-\bar{X})}^{2}}\sqrt{\displaystyle {\sum }_{i=1}^{n}{({Y}_{i}-\bar{Y})}^{2}}},\end{eqnarray}$$ ()

Julius et al.^[33,34] obtained three atomic models of the TRPV1 channel in different states (closed, half-open and open state) using high-resolution single-particle cryo-EM. Analysis of these structures reveals that there are two constrictions in the ion permeation pathway, which locate in the selectivity filter (denoted as the upper gate) and the lower gate region (Figs. 1(a) and 1(b)).^[33,50–52] The upper gate is formed by G643 at selectivity filter and the lower gate is formed by I679 at S6 helix.

To investigate the conformational changes of TRPV1 at different ambient temperatures that lead to the opening of the channel, we performed MD simulations at different temperatures based on the closed-state atomic model of TRPV1 (see Methods). The root mean square deviations (RMSDs) of the simulations (100 ns) are shown in Fig. 1(c). It is predictable that the RMSD becomes larger with the increment of temperature, especially in 410 K. Analysis using the HOLE program^[53] shows that the open-state conformations of TRPV1 are captured in the MD simulations at all temperatures calculated (Fig. 1(d) depicts a representative structure of open-state TRPV1 channel obtained from MD simulations). Inspection of the MD simulation trajectories shows that the open events of the upper and lower gate are independent of each other.

To show the state of the upper gate with time evolution, the time course of the channel open events of the upper gate at G643 is exported with 1 denoting the open state and 0 denoting the closed state (Fig. 2 and Fig. S1, see Methods). From Fig. 2, it is obvious that, after the equilibrium run, the constriction of the upper gate opens and maintains in the open state in most of the simulation time. This characteristic even appears at 310 K, which is considered to be unable to activate the channel. This result indicates that the opening of the upper gate is temperature-independent and the open state is a favored conformation for G643.

It is worth noting that, at temperature of 410 K, the upper gate is in the closed state in most of the simulation time (83%). Inspection of the simulation trajectory shows that a sodium ion occupies the pore at the G643 position and interacts with all the four G643s’ oxygen atoms of the tetramer (Fig. 3). The interactions pull the four G643s together and close the upper gate. We speculate that this phenomenon of ion blocking the pore is an occasional event because it is absent in other simulations in which cations only interact with one or two G643 residues. This result is consistent with that reported in the literature,^[52] which revealed that the carbonyl oxygen atom of G643 is the cation-binding region of the inner entrance. The interactions of cations with carbonyl oxygen atom of G643 may be important to the ion selectivity by attracting cations and repulsing anions. Inspection of the trajectory shows that after occupying this location ∼ 40 ns (from 10 ns to 50 ns), the sodium ion detaches from G643 and another sodium ion occupies this position (Fig. 3).

The cryo-EM structures illustrate that the residue M644 adjacent to G643 in the pore region has two different conformations in channel’s closed (PDB ID: 3J5P) and open state (PDB ID: 3J5Q). Both the states of M644 cannot form constriction to ion permeation. Statistical analysis of the simulation results reveals that the side chain of M644 prefers to locate in a hydrophobic pocket formed by I642, the hydrophobic portion of side chain of K639, and F640, Y671, V667, L647 of adjacent chain (Fig. 4(a)). This conformation is consistent with that shown in the cryo-EM structure (3J5P) considered in the open state. We chose the distance between the methyl group of M644 and C_α of Y671 as the marker to calculate the probability distribution of M644’s side chain (Fig. 4(b)). From the distributions in different temperatures, simulations in 310 K and 340 K clearly show the two states of M644. With higher temperatures, this two-state distribution becomes not obvious. The side chains of M644 s in the simulation at 360 K temperature stay outside of the hydrophobic pocket in all the simulation time. The initial conformation of all the simulations has the side chain of M644 being outside the hydrophobic pocket, therefore the entrance of M644 to the pocket is a spontaneous event in our simulations and seems to be temperature-independent.

In the simulations with temperatures 310 K and 360 K, M644 appears to have another conformation with the distance between M644 and Y671 in the range of 16–18 Å. In this state, the side chain of M644 bends to and blocks the pore. Inspection of the corresponding conformation reveals that bending of only one of four monomers’ long side chain of M644 can totally block the pore (Fig. S2). It appears that blocking of the pore by M644 is occasional and temperature-independent. In these three states of M644, the states with side chain of M644 buried in the hydrophobic pocket and bent to the pore region are relatively stable. Within the limited simulation time (100 ns), the conformational change between these two states is not observed.

From the above two sections, although G643 and M644 may occasionally block the pore, the blocking events are temperature-independent. In most of the simulation times, G643 is in the open state and cannot form a gate to restrict the permeation of ions. Considering that this position can form strong interactions with positive charges, maybe it is the origin of the permeation of cations for this kind of channels. M644 tends to be stably located in the hydrophobic pocket. In this state, it cannot form a gate to the ion permeation. However, the existence of a third state with side chain of M644 bending to the pore reveals that M644 can also influence the open probability of the channel.

The constriction of the pore at the “lower gate” portion is formed by I679. We also calculated the time course of open events at I679 (see Methods, Fig. 5 and Fig. S3). From Fig. 5, the open events of the lower gate at 310 K exist but rare. This is consistent with the experimental results that TRPV1 can occasionally open below the activation temperature. The amounts of open events increase with the increasing temperatures. This result indicates that the open of the lower gate is temperature-dependent. Figure 6 quantitatively shows the relationship between open probability of the lower gate and the temperatures. Calculations of the open probability are based on all the simulation trajectories (one 100 ns trajectory and two 30 ns trajectories). It is worthy to note that even up to 410 K, the open probability of the lower gate is 39.4% and cannot maintain in the open state (Fig. 5).

We investigate the mechanism underlying the temperature-dependent opening of the lower gate. Because the thermal effect is a factor impact the whole system, including the protein tetramer, the lipids and the surrounding water and ion molecules and small temperature difference (less than 10 °C) is sufficient to activate the protein, TRPV1 should not need very large conformational changes to open the pore. This can be seen from the comparison of the experimental cryo-EM structures in the apo and agent-binding state that opening of both upper and lower gate only needs the gating residue having a displacement of only 2–4 Å.^[33]

Opening of the lower gate needs I679 having an outward movement. Comparison of the lower gate’s closed and open state conformations obtained from the simulation trajectories shows that the major conformational change occurs at the N-terminal half (the lower half, from L577 to G590) of the S5 helix and the S4–S5 linker portion (from Q560 to R575) which moves toward the S4 helix when the lower gate opens (Fig. 7(a)). Inspection of the trajectories shows that the movement of S5 helix originates from formation of two hydrogen bonds (HBs) between E570 and S510/S512 (Fig. 7(d)). This result is consistent with the results of Wen et al.^[37] that formation of these two HBs is temperature-dependent.

Formations of these two HBs can dramatically torsion the S5 helix and the S4–S5 linker, which is unfavorable to the conformational energy. Therefore, both the HBs between E570 and S510/S512 are unstable and S5 and S4–S5 linker can return to the initial state and break the two HBs. In this way, the HBs between E570 and S510/S512 are in the dynamical switching between forming and breaking state. The dynamical switching between two states of the two HBs causes the S5 helix and S4–S5 linker have an outward oscillation. Principle component analysis (PCA) performed on S5 helix conformed that this oscillation movement dominates the movement patterns of S5 helix at all temperatures (Fig. 7(c) and Movies S1, S2). This result indicates that the outward oscillation of S5 is an inherent characteristic of TRPV1, which is consistent with the experimental result that TRPV1 channel can occasionally open below the threshold temperature.

It is proper to speculate that the oscillation of S5 helix should be influenced by temperature since temperature affects both the thermal movement of S5 and the formation probabilities of HBs between E570 and S510/S512. Since the oscillation movement is caused by the dynamical forming and breaking of the two hydrogen bonds, we calculate the formation rates of both the HBs (see Methods and Fig. 6). The formation rates of both the HBs show temperature-dependent character except for HB between E570 and S512 at 410 K. Temperature of 410 K is far beyond the endurable temperature of protein and can cause large thermal movement of protein’s secondary structure which replaces the effect of formation of the HBs in this instance.

We next investigated the relationship between the movement of S5 helix and the opening of the lower gate. Inspection of the structure shows that S5 and S6 helixes form tight hydrophobic contact within the same monomer (Fig. 7(b)). This hydrophobic contact ensures that outward movement of S5 can transmit to the S6 helix and cause same outward movement of S6 helix. To quantitatively show this correlation, we calculated the Pearson Correlation Coefficient of movements of the S5 and S6 helixes, which are represented by the outward motion of D576 and I679 respectively (see Methods and the locations of these two residues are shown in Fig. 7(b)). The calculated coefficients are 0.96, 0.92, 0.60 and 0.47 at temperature 310 K, 340 K, 360 K and 410 K, respectively. This result shows that the correlations between movement of S5 and S6 helixes are strong at 310–360 K, and that becomes relatively weak at 410 K. The simulation trajectory confirmed that, at temperature 410 K, the lower portion of the interface between S5 and S6 helixex cannot maintain the hydrophobic contact and water molecules entered this region due to strong thermal movements.

The above analysis shows a mechanical pathway that the temperature-dependent formation of two HBs between N-terminal portion of S5 and S6 helixes induces outward oscillation of S5 helix. Through the hydrophobic contact between S5 and S6 helixes, outward movement of S5 helix can transmit to S6 helix and open the lower gate. It is worth noting that the site of the formations of these two HBs is similar to the binding site of ligand (e.g., capsaicin and RTX) to TRPV1, in which the key residue is considered to be Y511.^{[26,33,39,54,55]} Our results indicate that the activation of TRPV1 by thermal and part of chemical ligand may originate from the same position which locates in the S3 helix.

Though our simulation results show a temperature-dependent opening of the channel pore, the exact open probability is not consistent with the experimental data (in which temperature of 340 K should fully open the channel). This is may be caused by the truncated structure used in this work. We added the missing residues in the pore turret part, but still lack of the C-terminal part of the TRPV1 channel. This part is considered to be the binding site of CaM and PIP₂, and can modulate the temperature sensitivity of TRPV1.^[56,57] We hope that the further work with the whole protein structure can show us the consistent data from theoretical calculations and experiments.

How heat activates TRPV1 channel is an important aspect of the gating mechanism of TRPV1 or even other heat-activation TRPV subfamily. Using molecular dynamics simulations, we investigate the channel-open process of closed-state TRPV1 at different temperatures (310 K, 340 K, 360 K and 410 K). The state of the proposed upper gate (G643) is temperature-independent and prefers to maintain in the open state. This site has the function to selective cation to pass the channel and is a main cation-binding site. M644 can occasionally block the pore through bend its side chain to the pore region. Oscillation of S5 helix originates from temperature-dependent forming, and breaking of two hydrogen bonds (E570-S510/S512) can induce the corresponding movement of S6 due to the hydrophobic contact between them. Outward movement of S6 helix releases the constriction of the lower gate to ion permeation. Clarification of the temperature sensor and heat-activation pathway of TRPV1 channel can help further drug discovery and functional study of TRPV1 channel.