Molecular properties of Ca2+ transport through TRPV2 channel: a molecular dynamics simulations study

Abstract TRPV channels are a category of nonselective cation channels that are activated by heat and ligands and permeate monovalent and divalent ions. The mechanism of Ca2+ transfer through TRPV2 channel is not well known. Here, we investigated the reaction coordination and energy fluctuation of Ca2+ transition in TRPV2 channel by steered molecular dynamics (SMD) simulations and potential of mean force (PMF) calculation. Results showed that electrostatic interactions between Ca2+ and residues of the first and second gates had main roles in ions transfer through the channel. Also, we recognized important amino acids in this path. Moreover, results indicated that enter and exit of calcium ions need to overcome barrier energies in the first and second gates. Communicated by Ramaswamy H. Sarma


Introduction
Transient receptor potential (TRP) channels are a heterogeneous superfamily of nonselective cation channels classified into six groups, including TRPA, TRPC, TRPL, TRPML, TRPP, and TRPV. The channels are regulated by various physical and chemical stimuli and inhibitors, which involve various cell actions (Muller et al., 2019;Per alvarez-Mar ın et al., 2013;Ramsey et al., 2006). Many TRP channels trigger and regulate signaling pathways via monovalent and divalent cationic fluxes (e.g., Ca 2þ , Mg 2þ , and Na þ ). TRP channels generally are similar to voltage-gated potassium or sodium channels in structural core as transmembrane topology (Yuan, 2019;Zheng, 2013).
Calcium is an important second messenger that regulates different signaling pathways and events in the cell, like gene expression, differentiation, proliferation, cell death, sensing, fertilization, etc. Changes in intracellular calcium concentration through channels (calcium influx) are a fast and localized process near the cytosolic mouth of the channels. Indeed, some physiological processes depend on a uniform and/or localized distribution of calcium in cells (Bikle et al., 2012;Mulier et al., 2017;Varghese et al., 2019).
The first Transient Receptor Potential Vanilloid (TRPV) family member is TRPV1, known to be the target of vanilloid compounds as capsaicin. They contain six members based on the similarity of amino acid sequence that the first four actions (TRPV1-TRPV4) relate to change of heat (Jin et al., 2006;Yuan, 2019). In structure, they are embedded in the membrane as tetramers so that each monomer consists of six transmembrane domains (S1 to S6). S1-S4 domains form a package that contributes to the thermosensing of channels, while S5 and S6 segments form the first gate (selectivity filter). Likewise, the distal end of S6 helixes forms a hydrophobic region in the second gate. Indeed, ankyrin repeat domain (ARD) is one of the most common motifs that are very conserved in TRPV channels. ARD consists of 6 ankyrin repeats that contain two antiparallel a helices followed by a finger loop (Cao et al., 2013;Chugunov et al., 2016;Zubcevic et al., 2016).
Here we focus on the TRPV2 channels, one of the interesting members of the vanilloid superfamily with widely gene expression in both neuronal and non-neuronal tissues. TRPV2 channel is a tetramer of four similar monomers (A, B, C, D). Each monomer contains 770 amino acids. From each monomer, 4 amino acids form the first gate (Gly606, Met607, Gly608, and Glu609). The second gate comprises the distal end of S6 from each monomer (AA630, AA650). In contrast to TRPV1, TRPV2 is not activated through capsaicin, but some reported compounds comprise 2-Aminoethoxydiphenyl Borate, Cannabinoid class (Cannabidiol, Cannabigerol, and Cannabidivarin, etc.) and Uricosuric compounds which active this channel (Koch et al., 2012;Per alvarez-Mar ın et al., 2013). Also, previous studies showed that thermal stimuli near 52 C can activate it and cause monovalent and divalent ions to pass across the channels (PCa 2þ /PNa þ ¼ 2.9). However, the thermal sensitivity of TRPV2 remains unknown in vivo. Therefore, further investigations are required to determine the function details of this channel (Caterina et al., 1999;Per alvarez-Mar ın et al., 2013).
The full-length structure of rabbit TRPV2 has been determined at high resolution by Cryo-electron microscopy ($4.4 Å). Huynh et al. in this study showed that the upper and lower gates of apo-TRPV2 are wider than the closedform of TRPV1 (Huynh et al., 2016). It seems that this structure is in a partially open state (Pumroy et al., 2019). Our study aimed to determine how calcium ions pass through the TRPV2 channel and recognize the important residues in this path.

Materials and methods
Preparation of simulation-system 3D structure of the TRPV2, was taken from the Protein Data Bank database (PDB ID: 5hi9). The Membrane Builder in CHARMM-GUI software was used to imbed TRPV2 into mixed lipid bilayers (POPC 256; POPE 121; CHOL 148) (Chugunov et al., 2016). In this system, 34806 TIP3 water molecules and 0.15 mol/L of NaCl (physiological ions concentration) were added to model the experimental conditions, and the net charge of the system was zero. So, the system was prepared for energy minimization.

Classical molecular dynamic simulation
This work performed all-atom classical molecular dynamics simulations using the GROMACS 4.5.6 package and CHARMM36 force field parameters. After minimization, NVT ensemble was done for the desired temperature of 310 K with periodic boundary conditions and using Langevin algorithm. Then, NPT equilibration was done for coupling pressure by Parrinello-Rahman barostat around 1 bar. Finally, MD simulation was run with a time step of 2 fs for 50 ns.

Steered molecular dynamic simulations
Umbrella sampling is a free energy calculation method that computes the potential of mean force (PMF) in some processes such as the binding affinity of ligand-receptor and passing ions or molecules through channels, etc. PMF indicates energy changes for any point of reaction coordinate through a series of conformational windows (Xu et al., 2020;You et al., 2019). So, to obtain the reaction coordinate path, we used SMD simulation based on time-dependent external forces. For SMD simulation, Ca 2þ was initially placed at the entrance of TRPV2 channel (with a 1 nm distance from the entrance of geometry of channel axis). Then, a constant force of 1000 kJ mol À1 nm À2 was applied to the Ca 2þ ion so that the ion did not move during energy minimization and 1000 ps NPT equilibration. Thereafter, for pulling Ca 2þ ion in the Z direction, a constant force of 1000 kJ mol À1 nm À2 was applied, while no constant force was imposed on the Ca 2þ ion. So, Ca 2þ ion was pulled parallel to the channel axis at constant velocities of 0.001 nm/ps from the mouth to the end of the channel. Ca 2þ ion movement was assessed relative to the reference group that consisted of COM of 10 amino acids (480-490 of monomer A). The total simulation time for the SMD was about 6 ns for the system. To calculate the PMF, the umbrella sampling simulations employed the same reaction coordinate as the SMD simulations, the z distance along the channel axis, and the reaction coordinate space is divided into windows every 1 Å. As a result, a total of 60 windows ranging from 0 to 6 nm in the z direction were individually simulated. Next, an NPT equilibration in each window was run during 1 ns, and then for each extracted input file, md-run was done for 5 ns. Finally, we constructed the PMF based on these umbrella sampling simulation trajectories using the Weighted Histogram Analysis Method (WHAM).

Rmsd
After embedding TRPV2 in the membrane, molecular dynamics (MD) simulation was done for 50 ns. The root mean square deviations (RMSD) value is a suitable parameter for evaluating structure stability because it makes the structural comparison between the initial structure and the following structure during simulations (Kufareva & Abagyan, 2011). As shown in Figure 1a, the system reached equilibration with RMSD value around 0.55 nm, indicating overall stability during MD simulations. After that, the RMSD value of steered molecular dynamics (SMD) simulation was calculated during 6 ns stabilized around 0.4 nm. The results indicated that the channel was stable when calcium passed through it ( Figure  1b). RMSD of MD and SMD were calculated for the backbone of TRPV2 channel.

Ion conduction pathway
The molecular architecture and gates region of TRPV2 channel has been described in some previous studies. The length of TRPV2 ion-conduction pathway from the beginning first gate to the end of second gate is approximately 50 Å (Huynh et al., 2016;Jin et al., 2006;Zubcevic et al., 2016). As shown in Figure 2a and b, the pore dimensions of TRPV2 were calculated after 50 ns simulations. The results indicated that the radius of the first gate had been wider than the second gate and showed a big cavity (about a radius of 6.5 Å) between the first gate and the second gate. The length of the first gate was about 10 Å that was located in a range of 9 to19 Å and consisted of four residues (Gly606, Met607, Gly608, and Glu609). Also, the length of the second gate was about 12 Å that located in a range of 32 to 44 Å in the ion conduction pathway. The second gate contained the distal end of S6 that Met640 and Ile642 were identified as important residues (Zubcevic et al., 2016). Previous studies showed that the first and second gates are an open state with Van-der-Waals radius of about 3.5 Å and 3 Å, respectively (Huffer et

Radial distribution functions
Radial distribution function (RDF) is a function of distance from a specific point that describes the density of atoms surrounding the point in MD simulation (Brehm et al., 2020). RDF was calculated to investigate the important residues in the transport mechanism of calcium ion through the TRPV2 channel. The average RDF in each umbrella sampling trajectory was calculated in equilibrium condition, and then the RDF graph was plotted along the 6 ns SMD. As shown in Figure 3a and b, RDF showed a sharp peak in monomer C of the first gate and monomer A and D of S6 segments (AA630-AA650) that confirmed ion passed close to these monomers. Moreover, results showed that maximum g (r) in the first gate was significantly higher than S6 segments (about 100 and 80, respectively) because the monomers of the first gate had an angle to Z direction (Figure 3c and d). On the other hand, by placing the ion at a distance of around 1 nanometer from other atoms, they can interact by electrostatic and Van-der-Waals interactions. So, these results showed that all monomers in the first gate and monomers A and D in S6 segments could be considered effective regions for passing ion through the channel.

Interactions energy analysis
The interactions energy was calculated for each residue in the ion transport path to determine which residues were important for passing calcium ion through the TRPV2 channel. Non-bonded energies that describe interactions between different atoms arise due to Van-der-Waals and charged interactions (electrostatic) modeled using the Lennard-Jones potential and Coulomb's law, respectively. Non-bonded energies between two stationary charged particles are calculated with the below formula in equation (1).
Where A and B are Lennard-Jones parameters, r is the interatomic distance, qi,j are the atomic charges, and D (¼ 4er) is a dielectric function. Moreover, interaction energy can be attractive or repulsive and shown as negative and positive values, respectively (Knegtel & Grootenhuis, 1998;Monasse & Boussinot, 2014;Wang et al., 2020). Indeed, the net interactions energy between atoms is obtained from the sum of these energies. In the first gate, while Lennard-Jones energy was positive in monomers A, C, D but Coulomb energy was highly negative. These results indicated that the interactions energy between calcium and the first gate was attractive. In addition, in the presented residues in the calcium ion pathway (S6 segments; AA630 -AA650), Lennard-Jones and Coulomb energies were calculated so that these energies in the monomer D were positive-negative, respectively, and the net energy of interactions was negative. Comparing these results showed that the attraction force (determined through interactions energy) between calcium ion and the first gate was more than in other areas (Figure 4).

Hydration of water in TRPV2 channel
In many channels, the hydration of water molecules plays an essential role in the transfer of ions (Jard on- Valadez et al., 2018;Richards et al., 2012). So, water density, hydration profiles, interactions of water with calcium ion, and diffusion coefficient of water (D w ) were determined during SMD simulation. Average of water molecules around first gate, S6 segments, and calcium ion were calculated in each umbrella sampling trajectory and equilibrium condition. Depending on the time of each window, the average of water molecules was plotted in time during 6 ns SMD. Our results showed that water density in the first gate was lower than the created hole by S6 segments (Figure 5a and b). Also, water hydration profiles indicated that the number of waters in the first gate fluctuates during the simulation as well as the second gate. However, in the second gate, the number of water molecules decreases dramatically from 4000 to 4300 ps compared to other times (Figure 5c and d). As shown in Figure 6b, calcium ion is located in the second gate region in this term. Therefore, it seems that dehydration has occurred due to the passage of calcium ions. To investigate the suggestion, hydration profiles of calcium ion was examined during SMD. This result showed that the number of waters around calcium ion decreased in this term (Figure 5e).
Moreover, the mean square displacement (MSD) method was also used to measure D w inside the channel using the Einstein relation. Using this method, we investigated the displacement of water molecules in the Z direction relative to the reference group (Wong-Ekkabut & Karttunen, 2016). Results indicated that D w inside the channel was 4.7198 ± 0.5655 Â 10 À6 cm 2 /s during the SMD simulation. In addition, the D w inside the channel, it was calculated for the first gate and the S6 segments too. Results indicated that D w in first gate was more than S6 segment so that their values were about 3.2856 ± 0.2586 Â 10 À5 cm 2 /s and 2.8965 ± 0.3850 Â 10 À6 cm 2 /s, respectively. For this reason, the water density in the first gate is lower than the S6 segments.  Furthermore, the Lennard-Jones energies between calcium ion and water molecules were positive in the first gate, and its values were lower than S6 segments. Unlike the Lennard-Jones and Coulomb energies were negative in the first gate, and their values were higher than S6 segments (Supplementary Figure S2). Moreover, fluctuation of Coulomb energy during 4000 to 4300 ps can be caused by the change of radius in this region, so the radius of this region was found to decrease from 6.5 Å in the middle channel to 2.5 Å near the end of the channel (second gate region) ( Figure 5). Moreover, as mentioned earlier, calcium ion dehydration occurred this term.
Results in Figures 4 and 5, demonstrated that interactions of water molecules with calcium ion caused decreased interactions of Ca 2þ with residues in the ion pathway. So, the interaction of Ca 2þ with residues in the first gate was more than the rest of the ion path

Ion permeation in TRPV2 channel
Constant velocity pulling is a method in SMD simulation in which an atom, molecule, or compound is moved with a constant velocity by applying external forces. The force is calculated below the formula in equation (2).
where k is the spring constant, x and v are the pulling atom's coordinate and velocity, respectively, and t is time (Hazel & Gumbart, 2017;Liu et al., 2006). If we want to keep the velocity constant, the force must change during simulation because some interactions affect the force. On the other hand, the force indicates whether the interactions are attractive or repulsive in the direction of coordination (C elerse et al., 2019;Sieradzan & Jakubowski, 2017). As shown in Figure 2b, the first gate was located 1 nm away from the channel's mouth, and its length was about 1 nm. Also, Figure 6b, showed that the calcium ion had traveled this path about 1.7 ns. Moreover, the force increased to pass ion thought first gate because of steric effect and abundance of interaction between calcium ion and residues. When the ion passed through the first gate, the force slumped, and calcium ion traveled a distance of about 1 nm during a few times step because the radius in this region (middle of S6 segments) increased to 6.5 Å, and the interaction between calcium ion and its amino acid decreased. Furthermore, as mentioned earlier, effective residues in this path and second gate caused to fluctuate in the force (Figures 4 and 6). Generally, these results indicated that required forces for passing calcium ion through the first gate were more than the second gate, which their values were about 1000 and 400 KJ/mol/nm, respectively.

Potential of mean force (PMF)
PMF is an essential parameter in computational studies that reflects all interactions of the desired molecule with its environment as ion, amino acids, solvent, etc. Moreover, dehydration at the entrance and inside the channel strongly affects PMF. It should be noted that PMF is not an average of interactions and potentials but is equivalent to the Gibbs free energy (Doudou et al., 2009;Sippl et al., 1996).
A big energy jump occurs when calcium ions enter from the first to the second gate. (Figures 1 and 6). Therefore, as shown in Figure 7, there were two high barriers at the first and second gates. In the first gate, energies for two main steps of the reaction mechanism were about 41.84 and 62.76 KJ/mol. Furthermore, as seen in Supplementary Figure  S1, there are two negative Coulomb interaction energies between the first gate and calcium ion in two-time intervals so that monomer A is effective to pass calcium ion in the interval of time from 0 to 1400 ps and monomer D and C have essential roles in it in the interval of time from 1400 to 2500 ps. in the second gate, a leap of energy occurred from 3.2 nm to 3.8 nm with an energy value of about 117.15 KJ/ mol and another energy step was about 12.55 KJ/mol. The obtained results indicated that barrier energy in the second gate was higher than the first gate because, in addition to the interactions energy, ion dehydration occurred at times 4000 to 4300 ps in the second gate (Figures 5e and 7). Also, the radius of the second gate was less than the first gate that was an important parameter in the energy calculation. Ultimately, the high energy barriers at the gates support this suggestion that calcium ion cannot easily enter and/or leave the TRPV2 channel.

Conclusion
SMD simulation of calcium ion pass through TRPV2 channel was used to analyze the conduction mechanisms. Our results show that there are four significant energy barriers in the calcium ion passageway, two steps in the first gate and the other two in the second gate, due to the interaction between calcium ion and residues in the first and second gates, as well as the decrement of channel diameter in these areas. Indeed, the transfer of calcium ion needs to overcome these deterrent forces. The forces in this path are of electrostatic type. The energy barrier in the second gate is higher than the first gate because, in addition to attractive forces between calcium ion and the residues, ion dehydration occurs. Also, the radius of the second gate was less than the first gate as an effective parameter in energy calculation. Indeed, transport of hydrated calcium ion through the first gate and dehydrated form in a narrow pore of the second gate are essential parameters in the conductance of ion in TRPV2 channel. Moreover, each residue's role in the conductions is determined by calculating Lennard-Jones and Coulomb energies. Generally, this study shows that calcium ion cannot enter or exit the TRPV2 channel easily.