Introduction

Ion channels are integral membrane proteins that control the flow of ions across cellular membranes. They play crucial roles in numerous physiological processes, such as the transmission of neural signals and the initiation of muscle contractions, among other functions. These channels can be categorized on the basis of their ion selectivity, distinguishing them as potassium, sodium, calcium, or chloride channels, or by their gating mechanisms, which include voltage-gated, ligand-gated, and mechanically gated types1,2,3,4,5. Dysfunction of ion channels can lead to a variety of diseases, often referred to as channelopathies, including neurological disorders, cardiac arrhythmias, and muscular dystrophies1.

Recently, a series of pioneering studies by Chang et al. demonstrated that terahertz (THz, 1012 Hz) electromagnetic fields can modulate the conductance of various ion channels, thereby influencing neuronal signaling6,7,8. These findings are highly significant and may pave the way for new therapeutic approaches to treating channelopathies. For instance, an experimental study revealed that a 5.6 μm wavelength (53.53 THz) field increased potassium current and modified action potential waveforms7. Another study showed that 8.6 μm wavelength (34.88 THz) field improved auditory perception in guinea pigs9. Simulation studies have also indicated that 51.87 THz field can increase ion permeation of the potassium channel KcsA10,11, and 53.60 THz field had a similar result on a model K+ channel12. 48.2 THz field was found to increase the ion permeability of sodium channels13. 42.55 THz field was found to alter the free energy profile of Ca2+ in the CavAb calcium channel6. 36.75 THz field was found to reduce the ion permeability of Kv1.2 by decreasing the stability of the channel, while the 37.06 THz field had the opposite effect14. 15 THz field has been reported to increase the ion permeability by increasing the proportion of direct knock-on permeation15,16. A lower frequency of 0.1 THz field was also found to alter the energy profile and accelerate the permeation of K+ in KcsA17. Some other frequencies, such as 1 THz, were also examined18. These findings suggest that THz fields may have direct regulation effects on ion channels. Beside ion channels, frequency-specific THz stimulation was also found to enhance surface wettability and induce water flux through artificial channels by energy transfer to water19,20. In addition, terahertz waves have been explored for their ability to alter gene expression, detect tumors, and potentially reduce side effects of drugs21,22,23,24, further demonstrating the therapeutic potential of the THz field.

Although the effects of THz fields on ion channels have been gradually recognized, most of the previous studies focused on the effects of THz fields on certain chemical groups of ion channels. Considering that there are usually well-defined ion binding sites within the selectivity filter (SF) of ion channels, one would expect that these bound ions should have well-defined oscillation motions, and the impact of THz fields on such motions remains largely uninvestigated. Therefore, this study aims to identify previously unknown frequencies that influence ion channel conductance, with a particular focus on those that directly affect the bound ions at the SF of ion channels.

We examined the structure and function of two well-studied ion channels: the voltage-gated potassium channel (Kv) and the voltage-gated sodium channel (Nav). The Kv channel is a tetramer with each subunit featuring an N-terminal domain, six transmembrane helices (S1 ~ S6) and a C-terminal domain. The voltage-sensing domain (VSD) is formed by S1 ~ S4, while S5, S6, and the intervening pore loop constitute the pore region. The SF, which controls the selectivity to ions, locates within the pore loop. It has a characteristic sequence TVGYGD, and provides five binding sites for K+ ions, named S0 ~ S4. The lower part of S6 helices form the gate, which undergoes conformational changes between different states1,3,4,5,25,26,27,28. The structure of a chimeric protein, Kv1.2/2.1 (PDB: 2R9R), resolved in 2007, has been extensively studied and serves as a basis for our open-state Kv1.2 model29,30,31,32,33 (Fig. 1A).

Fig. 1: The simulation systems.
figure 1

(A) Kv1.2, (B) Nav1.5. Purple and yellow spheres represent K+ and Na+, respectively. Each left panel shows the structure of the selectivity filter with the bound ions, and each right panel shows the pore radius along the z-axis.

Full size image

The eukaryotic Nav channel comprises an α subunit with four similar domains (DI ~ DIV), each containing six transmembrane helices (S1 ~ S6), along with several auxiliary β subunits. The SF is located at the pore loop between S5 and S6, formed by DEKA and EEDD residues from four domains, and the gate is formed by the lower part of S6 helices1,2,34,35,36,37,38,39,40,41,42. Nav1.5 is a sodium channel expressed in cardiac muscle and encoded by SCN5A43,44, whose partially-open structure was solved in 2021 (PDB: 7FBS)41. A recent study by Choudhury et al. proposed that the channel fully opens when the S6 helices in domains DI, DIII, and DIV or all four domains form partial π-helices45. Therefore, the Nav1.5 structure model with four manually-introduced π helices was used for this study (Fig. 1B).

In this study, we performed molecular dynamics (MD) simulations to investigate the oscillation motions of bound ions in the SF of Kv1.2 and Nav1.5, as well as the effects of THz frequencies on the Kv1.2 and Nav1.5 channels. Our findings revealed that THz waves with specific frequencies, 1.4, 2.2, and 2.9 THz, enhanced the ion permeation of Kv1.2, and 2.5 THz field enhanced that of Nav1.5, by directly affecting the motion of K+ or Na+ bound at the SF of the channels. 48.6 THz field increased the ion permeation of Nav1.5 by affecting the carboxyl groups (-COO−) at the SF. Conversely, 10.8 THz field decreased ion permeability of Kv1.2 by altering ion distribution. These results further suggest that THz fields could potentially be harnessed for therapeutic applications, offering a novel avenue for manipulating ion channel function.

Results

Structural analysis

We first computed the pore radius of each ion channel using HOLE46 after equilibration to confirm that they adopt a stably open state (Fig. 1 right panels). The pore radius of Kv1.2 was 0.60 Å at SF and 2.90 Å at the gate. The pore radius of Nav1.5 was 2.07 Å at SF and 2.55 Å at the gate. Both channel models were in their open state, and the SF was the narrowest location in the path of ion permeation.

Frequency analysis

To identify potential frequencies that influence the motion of the ions within the channels, we performed a frequency analysis using Fast Fourier Transform (FFT) on the trajectories of ions and the atoms or chemical groups that form the ion binding sites within the selectivity filter. This analysis encompassed a frequency range from 0.1 to 100 THz. The ions remained at the binding sites throughout the simulation (Fig. S1). In the following paragraphs, the z direction is perpendicular to the membrane and parallel to the channel, and the x/y directions are parallel to the membrane.

Our findings identified distinct characteristic frequencies for both ions and the key atoms or chemical groups within the channels. In the case of the potassium channel Kv1.2, three potassium ions occupied the SF after equilibration, binding to S0, S2 and S3/S4 (designated as K1 to K3 respectively). They exhibited vibrational motions predominantly within the 1 to 6 THz range. Specifically, the peak frequencies for these ions in the z direction were 2.2, 2.9, and 1.4 THz, respectively, while their peak frequencies in the x or y directions were approximately 5.0 THz (Fig. 2A).

Fig. 2: Frequency spectra.
figure 2

Frequency spectra of (A) K+ ions within Kv1.2, with the inset showing the z-direction spectrum of the three individual ions. (B) Oxygen atoms at the selectivity filter (SF) of Kv1.2. (C) Na+ ion within Nav1.5, with the inset displaying its z-direction spectrum. (D) Carboxyl (-COO) groups at the SF of Nav1.5.

Full size image

Additionally, we analyzed the carbonyl oxygen atoms of residues 370 to 375 in Kv1.2, which form the ion binding sites. We also performed FFT on the bond, angles and dihedrals that the oxygen atoms involve, so as to identify the vibrational modes. The atoms showed a peak frequency of 51.0 THz, corresponding to the stretching vibration mode of the C=O bond, which has been examined in the previous studies10,11,12. Another peak was observed at 10.8 THz in z direction, which mainly resulted from the CA-C-O angle vibration. A third peak at around 3.9 THz in the x-y plane was the result of N-CA-C-O and CA+1-N+1-C-O dihedrals (+1 denotes the atoms from the next residue) (Fig. 2B, Figure S2).

For the sodium channel Nav1.5, the Na+ ion within the SF exhibited a broad frequency spectrum ranging from 1 to 12 THz. A prominent peak at 2.5 THz was associated with the ion’s motion in the x-z plane within the SF, while a smaller peak at 10.0 THz corresponded to motion in the y-z plane (Fig. 2C), showing a similar pattern as observed by Zhu et al.47.

We also computed the carboxyl groups of Glu and Asp residues that interact with Na+ ions at the SF (residues 373, 376, 901, 904, 1425, and 1716). These groups exhibited peak frequencies at 48.6 THz and 42.5 THz, both linked to the stretching mode of the C-O bond. The 48.6 THz mode predominantly occurred in the x-z plane, while the 42.5 THz mode was observed in the x-y plane. A third peak showed up at around 2.8 THz, close to the peak frequency of the Na+ ion (Fig. 2D).

The effect of THz electromagnetic field

We initially performed control experiments without applying THz fields to establish a baseline for ion permeation, measured by the channel conductance. The conductance was calculated by counting the number of ions passing through the channels over a specified time under certain transmembrane potentials.

For Kv1.2, the simulated conductance was 12.6 ± 3.2 picosiemens (pS), consistent with the experimental range of 14 to 18 pS25. Throughout these simulations, ions permeated stably through the channel, with an average of three ions occupying the SF.

The simulated conductance of Nav1.5 was 39.1 ± 12.3 pS, exceeding but close to the experimental range of 10 to 20 pS41,48. This discrepancy may be due to the use of a more open structure in our simulations. Despite the elevated conductance, ion permeation through Nav1.5 remained stable, with an average of one sodium ion occupying the SF during the simulation period.

Subsequently, we applied electromagnetic fields at the observed frequencies to each ion channel system with the focus on the z-direction frequencies, as they are more likely to affect ion permeation11. The THz field was applied in the z direction with an amplitude of 0.8 V/nm, and channel conductance was measured under each condition. To determine whether the THz field had a direct impact on the motion of ions, we calculated the difference in average kinetic energy (ΔEk) of the ions located at SF compared with ions outside the channel (i.e. in the solution). The trajectories were split into 50 ns intervals, and the kinetic energy was averaged within each interval to reduce fluctuations. Additionally, we computed the ion distribution and potential of mean force (PMF) along z axis to determine whether the observed effects were associated with changes in ion distribution and PMF.

1.4, 2.2 and 2.9 THz fields increased ion permeation of Kv1.2

For Kv1.2, applying electromagnetic fields at frequencies of 1.4, 2.2, and 2.9 THz significantly increased channel conductance to 31.4 ± 6.6, 27.2 ± 6.5, and 30.8 ± 9.3 pS, respectively (n = 6 per frequency, p < 0.01). The effects were similar across these frequencies probably due to their proximity, with the 1.4 THz field producing the highest conductance, representing a 2.5-fold increase (Fig. 3A, B). A significant increase in the kinetic energy of the ions was observed (1.16 ± 0.06, 2.24 ± 0.17, and 1.33 ± 0.17 kJ/mol respectively), indicating that the THz field indeed enhanced ion permeation by directly affecting the ions (Fig. 3C). We noticed that 2.2 THz field yielded the largest ΔEk, but had no advantage over 1.4 or 2.9 THz on ion permeation. Since 2.2 THz corresponds to the vibration in S0, which locates near the surface of the channel, one likely reason is that it is less obstructed and easier to absorb energy.

Fig. 3: The effect of THz field on ion permeation.
figure 3

(A, D) The effect of THz electromagnetic fields at different frequencies on the ion permeation of Kv1.2 (A) and Nav1.5 (D). (B, E) The time-resolved average number of ions permeating through Kv1.2 (B) and Nav1.5 (E). The curves were smoothed by taking the sliding average in 2 ns intervals. (C, F) ΔEk of ions at the SF of Kv1.2 (C) and Nav1.5 (F) when different THz fields were applied.

Full size image

10.8 THz field hindered ion permeation in Kv1.2 by altering ion distribution

In contrast, a 10.8 THz field decreased the conductance to 3.0 ± 1.0 pS (n = 3, p < 0.01), approximately one-quarter of the original value (Fig. 3A, B). No significant changes in ion distribution were observed with the 1.4, 2.2, or 2.9 THz fields. However, the application of a 10.8 THz field altered the ion distribution and PMF. The occupation state of K+ shifted to S1, S2/S3, and S4, as opposed to S0, S2 and S3/S4 without the THz field (Fig. 4A, B). By examining the Fourier coefficient of each atom at SF, we noticed that the O atom in residues 370, 371, and 373 were the main contributors to this vibrational mode. The O atom in residue 370 and 371 altered the energy barrier between S2/S3 and S3/S4, and the O atom in residue 373 changed the preferable binding site of K1 from S0 to S1. Therefore, it is likely that this alteration led to the decrease in ion permeation.

Fig. 4: The effect of THz field on ion distribution and PMF.
figure 4

(A, B) The distribution and PMF of K+ along the z-axis within Kv1.2. (C, D) The distribution and PMF of Na+ within Nav1.5.

Full size image

2.5 and 48.6 THz field increased ion permeation of Nav1.5

For Nav1.5, the application of a 2.5 THz electromagnetic field dramatically reduced the residence time of Na+ ions and increased channel conductance to 133.5 ± 61.8 pS (n = 6, p < 0.05), representing a 3.4-fold increase compared to control conditions (Fig. 3D, E). A ΔEk of 0.81 ± 0.19 kJ/mol was observed, indicating that the THz field directly affects the motion of Na+ (Fig. 3F). In contrast, applying a 10.0 THz field resulted in an unaffected conductance at 39.7 ± 12.6 pS, although it induced an even higher ΔEk of 2.18 ± 0.77 kJ/mol. This indicates that the kinetic energy was not the sole reason for the increase in ion permeation.

The 48.6 THz field also enhanced ion permeation, but this effect was preceded by a 140 ns low-permeability delay period. Given that a brief 140 ns is significantly shorter than most physiological processes, this was considered as a transient effect. Therefore, we extended the simulation to 600 ns, and analyzed the trajectory after 140 ns when ion permeation became stable (Figure S3). The conductance increased to 93.0 ± 22.6 pS (n = 4, p < 0.05). The Na+ ions at the SF exhibited a ΔEk of 1.14 ± 0.30 kJ/mol, which affirms that the motion of carboxyl groups can absorb and transfer the energy to Na+, as proposed by Zhao et al.13.

The free energy barrier was reduced when THz field enhanced ion permeation

Previous works discovered that THz field can alter the distribution and free energy profile of ions, which can possibly affect the ion permeability6,15. Here, we computed the ion distribution and PMF to identify whether the above THz field can enhance ion permeation by altering the energy profile.

When the 1.4 THz field was applied to the Kv1.2 channel, no significant differences in ion distribution were observed, indicating that the binding sites remained unchanged (Fig. 4A). Meanwhile, the energy barriers between the binding sites, defined as the difference between peaks and troughs in the PMF curve, was reduced (Fig. 4B). The largest energy barrier, through which an ion at S0 escapes to the solution, was reduced from 17.8 kJ/mol to 14.4 kJ/mol. Similar results were observed when 2.5 and 48.6 THz fields were applied to Nav1.5. The ion density in the SF was lower, and the energy barrier decreased from 10.9 kJ/mol to 6.9 and 8.2 kJ/mol respectively (Fig. 4C, D). As a result, the ions can overcome the energy barrier easier and permeate faster.

THz field reduced hydrogen bonds between protein and water

Furthermore, we calculated the number of hydrogen bonds between the protein and water with the focus of hydrogen bonds formed by SF. For both Kv1.2 and Nav1.5, the number of hydrogen bonds slightly decreased when THz fields were applied. In Kv1.2, the average number of hydrogen bonds formed between protein and water decreased from 298 to 290 in the presence of a 1.4 THz field, and decreased to 279 when a 10.8 THz field was applied (Figure S4 A). The number of hydrogen bonds between SF and water also decreased from 7.2 to 5.9 and 3.9, mainly affecting the hydroxyl group of residue Thr370 (Figure S4 B). For Nav1.5, the average number of hydrogen bonds between protein and water decreased from 538 to 504 and 490 when a 2.5 THz or 48.6 THz field was applied (Figure S4 C). The number of hydrogen bonds between SF and water also decreased from 23 to 19 and 18, respectively, mainly resulted from fewer hydrogen bonds at the carboxyl groups in SF (Figure S4 D). This agrees with previous findings that the THz field can destabilize hydrogen bonds formed by the channel and water15,16, although this effect does not necessarily concur with an increase in ion permeability.

Simulations under non-specific conditions

The effect of non-specific frequencies

To assess whether the observed effects were specific to the peak frequencies, we performed additional simulations using non-specific frequencies. For this purpose, we chose a 75.0 THz field, which was distant from the peak frequencies identified in our spectrum and thus less likely to cause interference.

No significant changes in ion permeation were observed when a 75.0 THz field was applied to each channel. The conductance values for Kv1.2 and Nav1.5 were 10.2 ± 3.7 pS and 36.9 ± 7.4 pS, respectively (Fig. 3A, D). These values were consistent with the control measurements, suggesting that the changes in ion permeation observed in the above sections were specific to those selected frequencies and not due to a general effect of the THz field. Additionally, no significant changes in ΔEk was observed, although the perturbation induced a larger fluctuation in ΔEk (Fig. 3C, F). Therefore, this control experiment reinforces the idea that the peak frequencies identified by FFT analyses have a specific effect on ion channels, rather than simply causing a general perturbation.

The effect of THz field in different directions

We further investigated the specificity of the THz field’s polarization direction on its effects by conducting simulations with the THz field applied in the x and y directions, in addition to the z direction.

For Kv1.2, a 1.4 THz field significantly affected ion permeation when polarized in the z direction but had negligible effects when applied in the x or y directions. Furthermore, although potassium ions showed a peak frequency of 5.0 THz in the x-y plane, applying an electric field at this frequency in the x or y directions did not significantly impact ion permeation. In contrast, applying the 5.0 THz field in the z direction increased ion conductance to 25.7 ± 2.1 pS (n = 3, p < 0.01), probably due to its proximity to the peak frequencies in the z direction (Fig. 5A). This increase was less pronounced than that observed with the peak frequencies, further supporting the frequency specificity of THz fields.

Fig. 5: The effect of THz field in different directions.
figure 5

The effect of THz field in different directions on the ion permeation of (A) Kv1.2 and (B) Nav1.5.

Full size image

Similarly, in the case of Nav1.5, applying a 2.5 THz field along the x or y direction did not result in significant changes in ion permeation. Interestingly, when a 48.6 THz field was applied in the x or y direction, the conductance notably decreased to 10.1 ± 3.1 pS and 17.5 ± 2.3 pS, respectively (n = 6 each, p < 0.05) (Fig. 5B). This decrease was accompanied by an increase in ion density and energy barrier at the SF (Fig. 4C, D), indicating a tighter binding of Na+ ions to the SF of the channel, which likely caused the reduced ion permeability. These observations indicate that the influence of the THz field on ion channels is indeed specific to its polarization direction.

The effect of THz field with different amplitudes

We further investigated the effect of THz field amplitude by conducting simulations with varying amplitudes. For Kv1.2, we applied a 1.4 THz electromagnetic field, while for Nav1.5, a 2.5 THz field was used. The amplitudes tested were 0.2, 0.4, 0.6 and 0.8 V/nm.

Both Kv1.2 and Nav1.5 showed a positive correlation between the field amplitude and ion permeation. The conductance of Kv1.2 increased at all three amplitudes, with a linear relationship between amplitude and ion permeation, evidenced by a correlation coefficient of 0.84 (Fig. 6A). Nav1.5 showed a similar but slightly different response. The conductance only increased insignificantly when the amplitude was 0.2 or 0.4 V/nm, while a greater enhancement in ion permeability occurred at the amplitude of 0.6 or 0.8 V/nm. This indicates a non-linear relationship between the THz field amplitude and the ion permeability, suggesting that a higher amplitude is required to significantly influence the ion permeability of Nav1.5 (Fig. 6D).

Fig. 6: The effect of THz field with different amplitudes.
figure 6

Relationship between THz field amplitude and ion permeation, ΔEk, and residence time of Kv1.2 (A–C) and Nav1.5 (D–F).

Full size image

ΔEk also exhibited a positive correlation with the amplitude for both Kv1.2 and Nav1.5 (Fig. 6B, E). This shows that as the amplitude of the THz field increased, the vibrational motion of the ions was stimulated to a higher extent, thus increasing the kinetic energy. We also noticed that ΔEk was not proportional to increase in conductance, and a different nonlinear effect appeared between Kv1.2 and Nav1.5, which could be due to differences in their structures.

The residence time of the ions within the channels adopted a negative correlation with the THz field amplitude, further confirming that the THz field influences ion motion and accelerates permeation (Fig. 6C, F). The residence time of K+ in Kv1.2 reduced from 63 ± 10 ns to 29  ± 6 ns, roughly halved when the amplitude increased to 0.8 V/nm. The residence time of Na+ in Nav1.5 experienced a greater decrease from 3.6 ± 0.4 ns to 0.7 ± 0.1 ns, about one fifth of the original value.

The ion density and PMF also exhibited a gradual change (Figure S5). In Kv1.2, the PMF peaks at z = 1.7, 7.2, 13.0, 16.2 Å decreased, and in Nav1.5, the PMF peak at z = −1.2 Å decreased and the well at z = 3.5 Å increased. Both resulted in a lower energy barrier and therefore facilitated ion permeability.

Discussion

Despite numerous experimental and computational studies on the effects of THz fields on ion channels, it remained unclear at the molecular level whether the observed effects were specific or simply due to thermal influence. In our simulations, the system temperature was consistently maintained at 310 K with a high thermal coupling rate. No significant temperature changes were observed, with fluctuations of only 1.75 K in the Kv1.2 system and 1.35 K in the Nav1.5 system. Moreover, if the effects of the THz field were purely thermal, they would not exhibit directional specificity. However, our simulations demonstrated that the THz field effects were both frequency- and direction-specific. Therefore, the observed changes in conductance in our simulations can be attributed to the non-thermal effects of the THz electromagnetic fields.

Our study demonstrates that the influence of THz electromagnetic fields on ion channels is highly specific to both frequency and polarization direction. THz fields resonating at the frequencies of potassium ions within Kv1.2 (1.4, 2.2, and 2.9 THz) and sodium ions within Nav1.5 (2.5 THz) directly enhance ion permeation by accelerating ion motion along the permeating direction (z-axis). In contrast, a 10.8 THz field, resonating with the carbonyl oxygen atoms of the selectivity filter, alters ion distribution within Kv1.2, resulting in decreased ion permeability. These findings suggest that ion permeation through ion channels can be modulated by THz electromagnetic fields through multiple mechanisms, underscoring a multifaceted approach to controlling molecular behavior. While most of the previous research focused on frequencies associated with the protein motions, which influence ion motion indirectly, our study also revealed frequencies that resonate directly with the permeating ions at the SF, yielding more pronounced effects. Interestingly, a very recent study on the open-state model of the KcsA channel discovered a characteristic peak of ion vibration parallel to the channel axis located at 2.4 THz49, which is very close to one of the peaks we found for Kv1.2.

Given the similar structure of other ion channels in the Kv and Nav families, we believe that our results can be generalized to other Kv and Nav channels. However, some key differences can be observed between the Kv and Nav channels. The K+ and Na+ ions exhibited close but distinct characteristic frequencies in the Kv and Nav channels, which resulted from the different structures of the SF and the different energy profiles of the ions. The two channels also differ in their responses to the amplitude of THz fields. The reason may be that Kv channels have more well-defined ion binding sites than Nav channels, which makes resonance vibration of K+ in the Kv channels more likely to occur than Na+ in the Nav channels.

Our study demonstrates that ion channel function can be modulated by specific THz fields, opening new avenues for applications in biomedical fields, particularly in the treatment of channelopathies. The non-invasive and non-destructive nature of THz waves, combined with their ability to be precisely directed and focused, presents a promising therapeutic modality. Moreover, this approach could be extended to identify specific frequencies that affect other biomolecules, offering a versatile tool for molecular manipulation.

Despite these promising findings, our study has some limitations. The simulated electric field amplitude of 0.8 V/nm is significantly higher than what might be safely applied in practical scenarios, and the simulation duration may not fully capture the effects over longer periods. Therefore, further experimental studies are essential to validate the impact of THz fields on ion channels and their effects on cellular behavior at the macroscopic level. These future investigations should focus on refining the intensity and duration of exposure to ensure both safety and efficacy in real-world applications.

Methods

Model systems

The Kv1.2 model was built with the open-state structure of a chimeric protein Kv1.2/2.1 (PDB: 2R9R). Constructed from rat Kv1.2 and Kv2.1, this protein has an identical pore region as Kv1.2, and produces a similar current to Kv1.229. The pore region (residues 309 ~ 417 of each subunit) and the K+ within were selected and embedded in 1-palmitoyl-2-oleoyl-glycero-3-phosphocholine (POPC) bilayer to mimic the cellular membrane environment. The system was solvated with an aqueous solution containing 150 mM KCl. The channel axis was orientationed along the z axis, and the membrane parallel to the x-y plane. The TIP3P model was used for water molecules. The simulation system was prepared using CHARMM-GUI50. The final system size was 7.66 × 7.66 × 8.26 nm3, containing 34 K+, 22 Cl−, 128 POPC and 8716 water molecules.

The Nav1.5 model was built likewise using the open-state structure with 4 Π-helices, constructed by Choudhury et al.45. The pore region (residues 236 ~ 430, 825 ~ 945, 1320 ~ 1478, 1643 ~ 1777) was selected and embedded in a POPC bilayer. The system was solvated in a solution containing 150 mM NaCl. The final system size was 8.33 × 8.33 × 10.84 nm3, containing 57 Na+, 42 Cl−, 147 POPC and 15907 water molecules.

MD simulations

We performed molecular dynamics (MD) simulations using GROMACS 2022.351, augmented with an in-house developed plugin that enabled the addition of a constant electric field in each spatial direction (x, y, or z). This capability allowed us to simulate ion permeation in the presence of both an oscillating electric field and a uniform electric field along the channel axis.

The simulations were conducted under various conditions (Table 1), each with 3 ~ 6 replicates. The control and specific conditions were simulated with six replicates, while non-specific conditions were repeated three times. A total of 24 μs simulation was performed in this study.

Table 1 Simulation systems, conditions, and replicates
Full size table

The CHARMM36m force field52,53 and periodic boundary condition (PBC) were employed for the simulation. Leap-frog algorithm was used with a 2 fs time step. Temperature was maintained at 310 K using the Nose-Hoover thermostat54, and pressure was kept constant at 1 atm using the C-rescale barostat55. Particle-Mesh Ewald (PME) was used for long-range electrostatic calculations56. Van der Waals forces were considered up to a cutoff of 1.2 nm, with a switching function applied from 1.0 nm to smooth the transition.

To maintain the open state of the channels, a position restraint of 100 kJ/(mol ⋅ nm2) was applied to the protein backbone, except for the selectivity filter region (residues 370 ~ 375 for Kv1.2, and residues 372 ~ 376, 900 ~ 904, 1420 ~ 1425, 1712 ~ 1716 for Nav1.5). A constant electric field E0 was applied along z direction to establish a +450 or −450 mV transmembrane potential, calculated as the product of E0 and the length of the box along z direction lz32

$$U={E}_{0}{l}_{{{\rm{z}}}}$$

The ratio of electric field to magnetic field in electromagnetic waves equals the light speed c. Given that the velocities of charged particles in the system are far lower than c, the impact of the magnetic component is negligible6. Consequently, only the electric component was considered in our simulations. The total electric field applied to the system was therefore

$${{\bf{E}}}(t)={{{\bf{E}}}}_{{{\bf{0}}}}+{{\bf{A}}}\cos (2\pi ft)$$

Where A and f denote the amplitude and frequency of the THz field.

Trajectory and frequency analysis

All trajectory analysis were performed utilizing MDAnalysis 2.7.0 package in Python 3.11. After the equilibration of each system, a 1-ns simulation was performed for the frequency analysis. The same simulation setup was adopted as above, except that no restraints or electric field were applied. The 1 ns trajectory was split into 20 intervals. For each interval, FFT was performed on the velocities of ions and atom groups using NumPy 1.25.2. The modulus of Fourier coefficients were averaged among all intervals to obtain their frequency spectra. The frequency spectra of bonds, angles and dihedrals are calculated likewise. The formula of FFT can be written as

$$A(k)={\Sigma }_{t}\,f(t){{{\rm{e}}}}^{-2\pi {{\rm{i}}}kt}$$

Where f (t) can be the velocities in each dimension, bond lengths, angles or dihedral angles.

The modulus of the Fourier coefficients in x, y, and z dimensions indicates the vibrational mode for a certain frequency. For example, if the coefficient in z direction is much greater than in x, y directions, this vibrational mode mainly occurs along z axis.

Calculation of conductance

The conductance of a channel was calculated from the number of ions passing through in the simulation time under the transmembrane potential45,57:

$$G=\frac{I}{U}=\frac{Nq}{Ut}$$

Where N is the number of ions passing through in time t, I and U denote the current and transmembrane potential, and q is the charge of each ion.

Calculation of PMF

The PMF was calculated from the logarithm of the ion density15,32,47,57:

$$E(z)=-kT\ln \rho (z)$$

where k and T are the Boltzmann constant and the temperature, and ρ(z) denotes the ion density inside the channel along the z-axis. Please note that, as our simulations were conducted under non-equilibrium conditions, the PMF should be considered as Pseudo-PMF.