Theoretical investigation on the effect of terahertz wave on Ca2+ transport in the calcium channel.

The question of whether terahertz (THz) waves can interact with ions in channels of nerve cells and cause a further reaction has attracted much attention. To answer this question, we investigate the spontaneous radiation generated by Ca2+ moving in calcium channels and the effect of THz radiation on the transport of Ca2+ by solving the mathematical physical model through Brownian dynamics (BD) simulations. It is obtained that the moving Ca2+ in a calcium channel can generate electromagnetic radiation, the corresponding spectrum of which is concentrated in the THz range. Meanwhile, both the ion number in the channel and the background temperature are proved to have significant effects on the spontaneous emission spectra. The studies also show that external THz radiation can accelerate Ca2+ transport through the ion channel. These results are expected to provide a theoretical basis for the future treatment of THz waves in the neurological field.

Introduction

Ca2+ and Calcium channels play important roles in many physiological and biological functions, from heart contraction to gene transcription, and from signaling in the nervous system to triggering muscle contractions. Therefore, it is important to understand the precise feedback mechanisms that allow calcium ions to enter cells.

It has been reported that in addition to traditional electrical stimulation, light stimulation (Kato et al., 2012), sound stimulation (Ibsen et al., 2015), pain stimulation (Caterina et al., 1997), and pressure stimulation (Ranade et al., 2015) can promote ion transport across the membrane. However, there are few reports on ion transport stimulated by high-frequency electromagnetic stimulation. THz wave, as a kind of electromagnetic wave with low electronic energy level and high temporal and spatial resolution, has been widely used in electromagnetic biology research at the cellular and even molecular level in recent years. The biophotons produced by various intracellular oxidative metabolism processes happen to be in the THz band (Miller and Dumas, 2010; Wetzel and LeVine, 1999), proving that organisms have a wealth of THz information. Thus, THz waves have various potential applications in biological physics: they can be used to analyze the fingerprint spectrum of biological materials, stimulate proton transfer in base-pair hydrogen bonds (Wang et al., 2020), and accelerate DNA molecular unwinding (Wu et al., 2020).

THz waves have also been proved to have great potential for application in the field of neural research. It was claimed that through experiments and simulations verification, nerve fibers can be regarded as dielectric waveguides transmitting THz and mid-infrared waves (Kumar et al., 2016; Liu et al., 2019; Zangari et al., 2018). Studies showed that THz waves can regulate neural activity, and can enhance the permeability of the Voltage-Gated Calcium Channel (VGCC) by influencing the corresponding chemical bonds of —COO- or —C=O(Li et al., 2021).

Based on available scientific reports, we need to think about whether high-frequency electromagnetic radiation will be generated during the oscillation and transport of ions in the channels, and whether the spontaneous radiation will be affected by external stimulus? Therefore, this paper uses Calcium ion channels as a research object to investigate the spontaneous radiation of ions moving in channels. The schematic diagram of the basic model is shown in Figure 1. The main contents of this paper are as follows: (1) based on previous studies on BD simulation of the ion channel(Corry et al., 2001; Hoyles et al., 1998; Li et al., 1998), the three-dimensional potential energy distribution of the calcium channel is obtained in this paper to study the electromagnetic radiation generated during the transport of calcium ions across the membrane, aiming to explain the conjecture proposed in the ref (Liu, 2018), which is “the existence of high-frequency electromagnetic signals from THz to infrared in the biological nervous system; (2) we investigate the influence of temperature and the number of Ca2+ on the THz radiation generated by the moving Ca2+ in the channel; (3) we study the regulation of the external THz field on the movement of single ion and multiple ions, aiming to reveal the preliminary theoretical demonstration of the influence of the THz field on the production and conduction of neural action potentials. Studying the effect of THz waves on the ion penetration in Ca2+ channels can provide an approach not only for correcting an abnormal Ca2+ conduction process for normal cells but also for inducing rapid apoptosis of unwanted cells like malignant tumors by increasing Ca2+ flow to overload calcium concentration(Orrenius et al., 2003).

Figure 1

Schematic diagram of THz waves regulating calcium ion transport across membranes and ions producing THz radiation

Results and discussion

Geometric properties of calcium ion channel

It is well known that the Ca2+ channel mainly consists of a central water channel, which narrows from the central water area to both sides and connects the inner and outer cell sap of the cell membrane (Corry et al., 2001; McCleskey and Almers, 1985). The highly charged glutamate residues are distributed on the channel proteins connecting to a segment of the membrane, forming a negatively charged ring, which acts as a selective filter for ion transmembrane transports.

The calcium channel model used in this paper is constructed based on the model in the Ref (Corry et al., 2001). The two-dimensional (2D) model of the Ca2+ channel is shown in Figure 2, where the diamonds represent glutamate residues, and the triangles represent the dipoles. The three-dimensional (3D) model can be obtained by rotating the 2D model 180° along the axis r = 0, so that there are four dipoles and glutamate residues distributed on either end of the 3D channel, respectively. The charge of each dipole is 0.6 × 10−19C, with the length of 5 Å at z = −17.5 Å. These dipoles are used to represent the charged side chains thought to form a ring around the entrance of the constricted region and their nearby counter-charges. It is noteworthy that the glutamate residues are distributed asymmetrically on the boundary of the neck area of the channel (Bahinski et al., 1997; Ellinor et al., 1995), which is used to simulate selective filtration. The four charges are placed at z = 10.5 Å, 11.83 Å, 13.17 Å, and 14.5 Å with adjacent angular distances of 90°. The glutamate residues are with the same charge of −1.3 × 10−19C for each and are placed on the inner wall of the channel with a distance of 1 Å from the boundary. The values of the charges in the model are obtained from Ref (Corry et al., 2001).

Figure 2

The physical model of calcium channel

(A) The 2-D model of the calcium channel, where ▲ represents the dipoles and ◆ represent the glutamic charges.

(B) The 3-D model of the calcium channel is obtained by the rotation of figure (A) by 180° along the axis.

The ion-channel interaction includes two components: a repulsive force because of the induced charges on the protein boundary and an attractive force because of the electrostatic interaction of the ion with charge residues and dipoles on the channel wall. The induced charges on the protein boundary are related to the ions' positions in the channel. To describe the relationship between the induced charge and the ion, we located the Ca2+ at an initial position (0, 0, -4 nm) firstly and calculated its induced charge on the protein wall. Then move the ions forward in steps of 1Å and calculate the induced charge on the wall of the channel protein at the next position. Because solving the Poisson equation at every time step will greatly increase computation, the look-up table method (Hoyles et al., 1998) is used in this paper to store the pre-calculated electric fields and potential values of ions at all positions in the channel.

The relationship between ions with the induced charge on the water-protein boundary is shown in Figure 3, and the calculation method is shown in the method details section. The induced surface charge increases as the Ca2+ comes close to the channel neck area. Meanwhile, the induced surface charge shows the same electrical properties as Ca2+, which will give the Ca2+ a repulsive force when it gets close to the boundary.

Figure 3

Polarization charge density on the channel wall

The x axis represents the position of the ions along the axis, the y axis represents the axial distance. The induced surface charge increases as the Ca2+ comes close to the channel neck area. Meanwhile, the induced surface charge shows the same electrical properties as Ca2+, which will give the Ca2+ a repulsive force when it gets close to the boundary.

Potential energy distribution

The force on ions in the channel mainly includes three parts: the repulsive force from the induced charge at the channel boundary, the attraction force from the negatively charged groups of the channel, and the Coulomb force among the ions. For single Ca2+, the potential energy profile is constructed by calculating the movement of ions from infinity to a fixed position. Repeating this process can obtain the channel potential energy diagram. The dielectric constant of water inside the VGCC channel is different from free water (ε = 78). Especially in the neck areas, because of the negatively charged groups, the dielectric constant of water is significantly lower than that of the vestibule and free water. As a compromise, the same dielectric constant of water (ε = 60) is used in the pores and the reservoirs. The calculation method of potential energy is shown in the method details section. To verify the correctness of the model, the results are compared with it in ref (Corry et al., 2001) as shown in Figure 4.

Figure 4

Comparison of simulation results and reference results

The solid red line is the potential energy curve on the channel axis without considering the fixed charge on the boundary, whereas the dotted red line is the calculation result in the corresponding literature; the solid blue line is the curve of potential energy on the axis of the channel when the interaction of negatively charged groups is considered, and the dashed blue line is the calculation result in the corresponding literature.

Ca2+ entering the channel from outside of the membrane will be constrained by the attraction force of the negatively charged ring near the lowest point of potential energy and will meet a steeply rising potential barrier that is proportional to the square of the Ca2+ charge before the exit. Figure 5A shows the potential energy diagram of the 2D potential energy along the section y = 0. The lowest potential energy appears near the channel neck and shows increasing trends from the central axis to both sides.

Figure 5

Potential energy distribution in the calcium channel

(A) The 2-D potential energy diagram of the calcium channel.

(B) The potential energy curve of the channel axis with different dielectric constants of water.

Ca2+ transport across the membrane in the form of hydration, but Ca2+ and charged residues at the channel boundary will generate strong local electric fields, making the polarization degree of water inside the channel significantly different from that of free water. The degree of polarization of water can be characterized by the dielectric constant. As water plays an important role in the Ca2+ transport crossing the channel, we calculate the influence of the dielectric constant of channel water on the barrier height. It is worth noting that the height of the channel barrier varies with the dielectric constant of channel water. Figure 5B shows the potential energy curve on the central axis varies with the permittivity of water. When the permittivity of channel water was chosen to be as 60, 50, and 40, the height of the barrier increased accordingly.

The quantitative relationship between barrier height and the permittivity of channel water is shown in Figure 6. The results show that as the dielectric constant of the channel water decreases, the barrier height increases significantly. Here, we make a compromise and set the dielectric constant of the water inside the channel to 60 to complete the following calculation.

Figure 6

The relation between the height of the channel barrier and the dielectric constant of the channel water

As the dielectric constant of the water inside the channel decreases, the barrier to overcome for Ca2+ to pass through the channel increases significantly.

Ion motion and spontaneous radiation analysis

In the simulation to obtain the dynamic equilibrium of Ca2+ inside of the channel, the initial configuration of the Ca2+ is randomly distributed in the simulation region by generating a stochastic number. Without an external electric field, the Ca2+ in the simulation area will be attracted by the negatively charged ring in the neck area of the channel and move to the lowest point of potential energy near 1.25Å and move randomly. By solving Equations (Equation 10), (Equation 11), the positions and velocities of the ion at different times can be obtained, as shown in Figure 7. The parameters used during the calculation progress are listed in Table 1.

Figure 7

The relationship between the position and velocity of the ion changes with time without external stimulus

(A) The axial position of single Ca2+. It will be constrained to do random movement near the lowest point of potential energy (z = 12.5Å) without external stimulus.

(B) The ion velocity varies with time in the directions of x, y, z.

Table 1

The relevant parameters of the Brownian dynamic method used to calculate Ca2+ motion

Parameters	Description	Value
mi	Ca2+ mass	6.6 × 10−26kg
Q	Ca2+ charge	3.2 × 10−19C
Di (Corry et al., 2001)	Ca2+ diffusion coefficient	0.79 × 10−9m2/s
KB	Boltzmann constant	1.38 × 10−23J/K
T	Temperature	310K
Ri	Ca2+ radius	0.99 Å
F0 (Stillinger and Rahman, 1974)	the force field parameter	2 × 10−10N
RB (Tansel et al., 2006)	Ca2+ hydration radius	2.6 Å
Rw (Stillinger and Rahman, 1974)	the radius of boundary wall atoms	1.4 Å
C0 (Corry et al., 2001; Guàrdia and Padró, 1996; Guàrdia et al., 1993)	Coefficient C0 for Ca2+-Ca2+ pair	0.8 KBT
C1 (Corry et al., 2001; Guàrdia and Padró, 1996; Guàrdia et al., 1993)	Coefficient C1 for Ca2+-Ca2+ pair	1.6 Å
C2 (Corry et al., 2001; Guàrdia and Padró, 1996; Guàrdia et al., 1993)	Coefficient C2 for Ca2+-Ca2+ pair	1.8 Å
C3 (Corry et al., 2001; Guàrdia and Padró, 1996; Guàrdia et al., 1993)	Coefficient C3 for Ca2+-Ca2+ pair	1 Å
cw (Corry et al., 2001; Guàrdia and Padró, 1996; Guàrdia et al., 1993)	Coefficient cw for Ca2+-Ca2+ pair	2.76 Å

Here, the spontaneous radiation spectrum of Ca2+ movement can be obtained by applying the Fourier transform of the electric field generated by calcium ion movement. The electric field emitted by ion motion can be obtained by Equation (1) (Jackson, 1999),where q is the charge of Ca2+, d is the distance between the measurement point and the ion, c is the speed of light, is the unit vector of the electric charge pointing in the direction of the observation point, and is the accelerated speed, which can be obtained by taking the derivative of Equation (11). The total radiation spectrum can be obtained from the radiation spectrum in the three directions. The results are shown in Figure 8. As can be seen, the spectrum of calcium ions moving in the channel is mainly concentrated around 16.9THz.

Figure 8

Radiation spectrum produced by Ca2+ in the calcium channel

(A) The radiation spectrum of a Ca2+ moving in a channel in the directions of x, y, z.

(B) The total radiation spectrum of a Ca2+ moving in the channel. The blue curve is obtained by Gaussian smoothing of the red curve, and the spectrum of Ca2+ in the channel is mainly concentrated around 16.9THz.

The process of ion transport is the continuous movement of ion current from one side of the membrane to the other. Based on the analysis of the motion of a single ion, we studied the changes in the spontaneous emission spectrum of the ion when there are two ions in the channel, the spectrum characteristics of which are shown in Figure 9. The results show that the center frequency of the radiation spectrum of the two ions remains the same, but the radiation intensity varies between the ions. Considering the superposition effect of energy, when multiple ions are transported through the membrane, a higher energy will be radiated outward, which can be measured by specific experimental methods.

Figure 9

In the calcium channel, the influence of the number of ions changes on the radiation spectrum of Ca2+

(A) Comparison of the radiation spectra of one Ca2+ and two Ca2+ in the channel.

(B) The red and blue curves are obtained by Gaussian smoothing of the data in (A).

Because the random thermal motion of ions is correlated with temperature, the effect of temperature on the ion motion is investigated and shown in Figure 10. We selected three temperature values (280K, 310K, and 340K) for this investigation. The dielectric constant of the channel water is assumed to be constant during the simulation. As the temperature increases, the radiation spectrum of Ca2+ shifts to higher frequencies. It can be understood that, with the increase in temperature, the thermal movement of ions intensifies, making the radiation spectrum shift to a higher frequency.

Figure 10

In the calcium channel, the influence of temperature changes on the radiation spectrum of Ca2+

(A) Spectral characteristics of Ca2+ at different temperatures.

(B) The blue, red, and yellow curves are obtained by Gaussian smoothing of the data in (A).

Applied THz electric field effect on single Ca2+ transport through the channel

The movement process of ions in calcium channels regulated by THz waves is studied by the BD method. The Schematic diagram is shown in Figure 11. In this section, we discussed the regulation of THz pulses during the transmembrane process of single ion and multi-ion systems and the changes of the spontaneous emission spectrum of ions under THz irradiation.

Figure 11

Schematic diagram of THz wave acting on calcium ion channels

The THz electric field waveform used in the simulation is a unipolar square wave, and the transmission direction is the r-axis, which is perpendicular to the axial direction of the ion channel; the electric field direction is the z axis, which is parallel to the axial direction of the ion channel.

Figure 12A shows the effect of unipolar THz pulse trains under different field amplitudes on the transport of single ions through the channel when the frequency is 1 THz. The position of the Ca2+ placed in the channel was calculated at each discrete time step of 2 fs for 2×106 time steps. Figure 12B shows the effect of unipolar THz pulse trains under different frequencies on the transport of a single ion through the channel. The electric field intensity used in this simulation is −6 × 108V/m. The results show that the time it takes for Ca2+ to cross the channel decreases with the increase of the amplitude. Meanwhile, the time it takes for Ca2+ to cross the channel decreases with the increase of the frequency of the THz wave.

Figure 12

Trajectories of ions irradiated by unipolar THz waves with different amplitudes and frequencies

(A) The influence of amplitude changes on Ca2+ motion when the frequency is 1THz.

(B) The influence of frequency changes on Ca2+ motion when the amplitude is −6 × 108V/m.

In addition, we compare the changes of ion radiation spectrum with or without an external THz field. Here, we focus on the influence of the THz electric field on ion transport across the membrane, while ignoring the influence of the magnetic field. The relevant analysis process can be referred to ref (Bo et al., 2020). We chose the unipolar THz pulse train as the excitation source. The frequency of the applied THz field is 1THz. The electric field waveforms in the time domain and frequency domain of the THz field are shown in Figure 13A, and the influence of THz waves on the spontaneous emission spectrum of Ca2+ is shown in Figure 13B. Under the irradiation of a unipolar THz field, the peak value of the ion's total radiation spectrum is more obvious at low frequency, and the peak position is the same as the peak position of the spectrum of the monopole THz electric field. Meanwhile, the effects of THz electric field irradiation on the radiation spectrum generated by ion oscillation at high frequencies are very weak.

Figure 13

Effect of THz electric field on the spontaneous emission spectrum of ions

(A) Time-domain and frequency-domain waveforms of the THz field used in the simulation.

(B) Radiation spectrum of ion motion under THz wave irradiation.

Applied THz electric field effect on multiple Ca2+ transport through the channel

In a practical calcium ion channel, millions of calcium ions pass through the channel every second. And the Ca2+ will be subjected to the repulsive forces of multiple ions to accelerate the penetration rate of ions. Figure 14 shows the comparison of the potential energy distribution when there are two ions in the channel with that of a single ion.

Figure 14

The potential energy on the axis of the channel

The blue line represents the potential energy curve of a single Ca2+. The red line represents the potential energy curve of one calcium ion under the influence of another ion in the channel neck. The coordinate of the calcium ion in the channel neck is (0, 0, 21.5Å).

Because of the Coulomb repulsion between ions and the effect of the induced charge of the second ion, the energy barrier is reduced. It indicates that the electric field used in the simulation of the multi-ion system is much lower than that of the single-ion system at physiological concentration. To obtain a better statistic, six Ca2+ are put on the right side of the channel. Without the external field, Ca2+ would be attracted to the glutamate groups to move to the vicinity of the channel neck. In this section, we mainly care about the movement of ions from one side of the channel to the other side. Thus, we assume that the movement of ions will not cause significant changes in the concentration of calcium ions on both sides of the membrane within the calculation of nanosecond magnitude. When an ion moves from the right to the left of the channel, a new Ca2+ will automatically fill the right side. The amplitude of the applied electric field is −3 × 108V/m, the variation of the number of ions across the channel with time in 8 ns is shown in Figure 15. The results show that as the frequency increases, the number of transmembrane ions increases significantly at the same time, which indicated that the high-frequency THz field was more conducive to regulating Ca2+ transmembrane penetration.

Figure 15

The variation of ion transmembrane number with time under the control of THz pulses of different frequencies

The amplitude of the applied electric field is −3 × 108V/m. The total calculation time is 8 ns, and the time step is 2fs.

Conclusions

In this work, we established a mathematical physical model of the Ca2+ channel by the Brownian dynamic method. The time domain and frequency domain characteristics of calcium ion movement in the channel have been analyzed. The results show that the Ca2+ transmembrane transport is a rapid signal change process on a short timescale, and its spontaneous radiation spectrum is mainly concentrated in the THz range. Meanwhile, the changes in temperature and the number of ions in the channel were proved to have significant effects on the spectral characteristics, as the frequency increases, the radiation spectrum of calcium ions shifted toward higher frequencies. In addition, under the external THz field irradiation, the oscillation spectrum of the irritated Ca2+ is significantly enhanced, which is expected to be applied in future nondestructive testing. It was also found from the simulation results that the transport process of Ca2+ across the ion channel is related to the frequency and amplitude of the THz wave. Under a certain amplitude, the ion permeability gradually increases as the frequency of the THz field increases; under a certain frequency, the ion permeability also increases as the amplitude increases. The finding in this paper provides a theoretical basis for the future treatment of THz waves in the neurological field.

Limitations of the study

The study is theoretical. The interaction between THz waves with ions in channels of nerve cells needs further experimental investigation.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to the lead contact, Yubin Gong (ybgong@uestc.edu.cn).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

Our study does not use experimental models typical in the life sciences.

Method details

BD calculation method

In BD (Gunsteren et al., 1981; Hoyles et al., 1996), the movement of a single ion through the membrane is determined by the Langevin equation of motion,where m is the ion’s mass; is the vector velocity of the ion at time t; γ is the damping coefficient which reflects the average frictional force under the action of water molecules around the ion; is the random forces; is the electric field force generated by ions and the applied THz field; represents the short-range forces, including collisions between ions and the repelling force of the channel boundary on the ion, which is used to maintain the stability of the system to avoid ions’ overlap and crossing of the water-protein interface.

In this simulation, we applied a rigid boundary on the inner surface of the channel to prevent ions from crossing the boundary (Corry et al., 2001). Therefore, the ion moving toward the channel’s inner surface will be repelled back. The short-range force between Ca2+ and the boundary can be expressed as (Guàrdia et al., 1991a, b),where F0 is the force field parameter, R is the radius of ion, R is the radius of boundary wall atoms, r is the distance between the ion and the central axis, and R is the distance between the channel wall and the central axis in the corresponding z coordinate.

For the ions moving in the channel, the forces between multiple ions are mainly described by the Coulomb force and the short-range force between ions. The short-range force (Guàrdia et al., 1991a, b) between two ions with separation r can be expressed as,where R1and R2 are the radii of the two ions. The meanings and values of coefficients C0, C1, C2, C3, c are listed in Table 1, which are provided from Ref. (Corry et al., 2001; Guàrdia and Padró, 1996; Guàrdia et al., 1993).

the Langevin Equation (2) can multiply the integrating factor and be integrated from an initial time t to t，

In addition to random forces, the forces exerted on the ion are denoted by . Then, we can get an expression for the ion velocity

Taylor expansion of at t,

Substitute Equation (7) into (6), and then we can get,

By integrating (8) from an initial time t to t, the position of the ion at the next moment can be obtained. Here, the result of integrating the last term of Equation (8) is expressed in terms of a defined X variable,Where, τ=γ(t-t). The expression without velocity term can be obtained from Equations (8) and (9),X(Δt) is the random variable used in the BD algorithm, which obeys Gaussian distribution with a mean value of 0 and a variance of 2KTΔt/mγ=2DΔt. One can refer to Ref. (Gunsteren and Berendsen, 1982) for details. Equations (10) and (11) are the displacement and velocity equations of ions, which are used to calculate the displacement and velocity of Ca2+ in the channel.

The total electric field force can be written as the following form,where is the force of induced charge generated by the induction of ions at the water-protein interface; is the force on the ion i generated by other ions in the channel; is the electric field force generated by the applied THz electric field. The electric field waveform used in the simulation in this paper is a unipolar square waveform. Here, since the motion velocities of ions are much smaller than the light speed, the effect of the magnetic field on the motion of ions is not considered (Bo et al., 2020). The values of the parameters are listed in Table 1.

Poisson equation solution

The total electric field at the ion position is determined by the solution of the Poisson equation,where ρ0 represents the charge density of the ion, ρ represents the charge density of the negatively charged residues at the channel boundary and the channel dipole. ε0represents the vacuum dielectric constant, which is 8.85×10-12. ε1 and ε2 represent the permittivity of water in the channel and channel protein, and the values are 60 and 2, respectively. The total potential in the channel can be expressed as follows,where φ is the self-potential due to the ion and the induced charge on the channel wall, φ is the external potential due to the negatively charged groups, and the dipole on the channel wall, φ is the external potential due to the applied field.

To get the potential distribution of the channel, we replace this model with an equivalent model having a dielectric constant ε1 that produces the same electrical potential throughout space. Furthermore, the discontinuity boundary between channel water and channel protein under the electric field is replaced by polarization charge density ρ1. Thus, the potential of Equation (13) can be expressed by,

By iteratively solving Equations (15) and (16), we can get the potential distribution of the ion channel. Here, is the normal direction of the boundary. The criterion for convergence is that when the ratio of the polarization charge density ρ1 difference between the two iterations to the maximum charge density is less than 0.001%, we consider the result to be convergent.

Channel potential energy calculation method

The energy difference caused by neglecting the change of the dielectric constant of water is approximately given by the Born energy,Where R is the ionic hydration radius, ε is the dielectric constant of free water, the value is 78. ε1 is the dielectric constant of water in the channel, and the value is 60. Here, a smooth switching function (Chung et al., 1999) is used to avoid complications arising from sharp potential energy changes in BD simulations, which can implement this potential barrier,

The potential energy (Chung et al., 1999) can be calculated by:

Quantification and statistical analysis

The BD program used in the simulations is written in MATLAB software and executed in the server with the Intel(R) Core(TM) i7-8700k CPU and 32GB memory. The algorithms have been deposited at Zenodo (Guo et al., 2021a, 2021b). The main calculation process includes two steps. First, the three-dimensional space is discretized into 11 × 11×44 points. The potential distributions of the channel for ions at different positions are pre-calculated which are stored in the 4-D data table, which has been deposited at Zenodo (Guo et al., 2021a, 2021b). The CPU time of a server for completing the whole calculation of the three-dimensional potential is about 6 h. Then, the Langevin equation of motion is solved and the positions and velocities of the ions at different times are obtained. By referring to the data table established in the first step, the electric field forces of ions at different positions can be obtained. The time step used in the dynamics simulation is 2 fs, and the temperature is 310 K. With one ion in the reservoirs, the CPU time of a server for completing a single simulation period of 400ps (2×106 time steps in 2 fs) is about 20 min.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Calcium channel model	This paper	https://doi.org/10.5281/zenodo.5728726
Software and algorithms
MATLAB	Mathworks	https://www.mathworks.com/
Regulation of calcium ion transport by THz wave algorithm	This paper	https://doi.org/10.5281/zenodo.5720258