Molecular dynamics simulation method based on action of terahertz wave and ion channel

Abstract

The invention relates to the technical field of biophysics, in particular to a molecular dynamics simulation method based on the action of terahertz waves and ion channels, and the molecular dynamics simulation method deepens spectrum-structure correlation analysis and expands the application range of ion channel types by optimizing a field effect simulation process. Technical support is provided for the fields of terahertz wave biological effect research, neural signal regulation and control mechanism analysis, novel physical treatment strategy development and the like, the technical bottleneck problem of a traditional method in research of electromagnetic field-biomacromolecule dynamic interaction is effectively solved, and the method has important scientific significance and clinical application prospects.

Description

Technical Field

[0001] This invention relates to the field of biophysics, and in particular to a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels. Background Technology

[0002] Ion channels, as core proteins in cellular electrophysiological activities, participate in signal transduction and regulate various important cellular activities, making them a key issue in life sciences and biomedical engineering. Terahertz waves lie between microwave and infrared frequencies. The vibrations of the biomolecular skeleton, such as proteins, and the weak interactions between molecules in biological tissues are located within this frequency range. Furthermore, due to the low-energy nature of terahertz waves, their use in regulating biomolecules does not damage the original molecular structure. Therefore, terahertz modulation holds great potential and application value in future physiological function regulation and disease treatment.

[0003] Currently, this interdisciplinary research field still lacks a systematic research paradigm. This is mainly reflected in the fact that the academic community has not yet established a standardized molecular dynamics method specifically for simulating terahertz wave-ion channel interactions. From the selection and loading of physical field parameters and the construction of simulation systems to the control of the computational process, there are no established rules to follow, resulting in a shortage of related research and severely restricting the exploration of its microscopic mechanisms. At the same time, existing technologies also lack effective analytical and detection systems for the specific biophysical effects that electromagnetic fields may induce. Traditional trajectory analysis tools are unable to accurately quantify field effects from atomic motion data, such as gating structure conformation and ion transport energy, making it impossible to establish a clear causal logic from physical stimulation to functional response. Therefore, this field urgently needs a complete technical solution that integrates standardized simulation processes with targeted analytical methods. Summary of the Invention

[0004] The purpose of this invention is to provide a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels, aiming to solve the problem that existing technologies cannot quantitatively evaluate the influence of terahertz fields on key indicators such as ion channel conformation, pore size, free energy, and ion binding.

[0005] To achieve the above objectives, this invention provides a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels, comprising the following steps: The target system was modeled to obtain a three-dimensional biomimetic structure of ion channel-phospholipid bilayer-solution environment; Based on the multi-stage energy minimization and ensemble equilibrium conditions, dynamic simulations of the tactile conformation of the system are performed to obtain the equilibrium system. By combining Fourier transform and linear response theory with the obtained equilibrium system, the vibration spectrum of the desired binding site is obtained; Based on the characteristic frequencies derived from the spectrum, terahertz waves of the corresponding frequencies are loaded and molecular dynamics simulations are performed again. The simulation data are analyzed in multiple dimensions to obtain quantitative indicators of the interaction between terahertz waves and ion channels.

[0006] In the section "Modeling the target system to obtain the three-dimensional biomimetic structure of ion channel-phospholipid bilayer-solution environment", an electroneutrality model can be established using CHARMM-GUI according to the required ion channel type, phospholipid type and liquid environment.

[0007] The section on "dynamic simulation of the system's thixotropic conformation based on multi-stage energy minimization and ensemble equilibrium conditions to obtain the equilibrium system" includes the following steps: By performing energy minimization on the established system, an energy-stable conformation is obtained; Using the NVT ensemble, the energy-stable conformation of the complex was equilibrated at temperature to obtain a temperature-stable simulation system; Using the NVT ensemble, pressure equilibrium was achieved in the temperature stability simulation system to obtain the equilibrium simulation system.

[0008] In the section "By combining Fourier transform and linear response theory with the obtained equilibrium system, the vibration spectrum of the required binding site is obtained", the vibration spectrum calculation mentioned refers to establishing the autocorrelation function of the total dipole moment of the entire system and performing Fourier transform and linear response theory to obtain the vibration spectrum of the system. By combining VMD, the functional groups of the system corresponding to the spectral peak can be obtained in real time, and the characteristic vibration frequency corresponding to the ion channel binding site can be obtained.

[0009] Specifically, in the section "Based on the characteristic frequencies derived from the spectrum, a corresponding terahertz wave is applied to perform molecular dynamics simulations again, and the simulation data is analyzed in a multi-dimensional molecular manner to obtain quantitative indicators of the interaction between the terahertz wave and the ion channel," the quantitative indicators include: The root mean square deviation of the conformational trajectory data is calculated to obtain the structural stability index of the complex. Based on the average potential calculation, the change in potential energy before and after the application of terahertz waves is obtained; The interaction energy of conformation trajectory data was analyzed using an energy decomposition algorithm to obtain the influence of terahertz waves on the system energy. Based on aperture calculations, the influence of terahertz waves on the aperture size of ion channels was obtained; Dynamic structural analysis of conformational trajectory data based on secondary structure reveals the influence of terahertz waves on the structure.

[0010] This invention presents a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels. It constructs a composite system model of a transmembrane protein-phospholipid bilayer-solution environment, enhancing the realism of the biomimetic simulation environment. Regarding the analysis of interaction mechanisms, the spectrum-guided multidimensional analysis method developed in this invention comprehensively analyzes the interaction mechanism between terahertz waves and ion channels from multiple levels, including conformational stability, free energy changes, pore dynamics, and structural retention. In particular, the combination of vibrational spectrum identification and linear response theory can accurately pinpoint the resonance frequency and quantitatively analyze the regulatory effect of the terahertz field on the gating and transport properties of ion channels, providing a theoretical basis for the physical field regulation of biomolecular functions. In terms of practicality and adaptability, the integrated simulation scheme provided by this invention can be flexibly adapted to various ion channel systems, such as voltage-gated ion channels and ligand-gated ion channels. Through standardized processes, it achieves full-link automation from system construction and spectrum analysis to field effect simulation, lowering the technical threshold for cross-disciplinary research, improving the reproducibility of electromagnetic biology research, and demonstrating its potential value in neuromodulation and innovative therapy development. The molecular dynamics simulation method of this invention provides technical support for the study of terahertz wave biological effects, the analysis of neural signal regulation mechanisms, and the development of novel physical therapy strategies by optimizing the field effect simulation process, deepening the spectrum-structure correlation analysis, and expanding the applicable scope of ion channel types. It effectively solves the technical bottleneck problem of traditional methods in studying the dynamic interaction between electromagnetic fields and biological macromolecules, and has important scientific significance and clinical application prospects. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0012] Figure 1 A schematic diagram of the established protein-phospholipid bilayer-solution environment system model.

[0013] Figure 2 This is a schematic diagram showing the vibrational frequencies of the protein selective filter structure calculated based on the established model.

[0014] Figure 3 The change in average potential (PMF) before and after loading terahertz waves.

[0015] Figure 4 The change in binding energy before and after loading terahertz waves.

[0016] Figure 5 The aperture change before and after loading terahertz waves.

[0017] Figure 6 The changes in the secondary structure before and after loading terahertz waves.

[0018] Figure 7 This is a flowchart of a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels provided by the present invention. Detailed Implementation

[0019] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0020] Please see Figures 1 to 7 This invention provides a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels, comprising the following steps: S1 models the target system to obtain a three-dimensional biomimetic structure of ion channel-phospholipid bilayer-solution environment; Electroneutrality models can be established using CHARMM-GUI, depending on the required ion channel type, phospholipid type, and liquid environment.

[0021] Specifically, based on the required ion channel type, phospholipid type, and solution environment, an initial complex system containing membrane proteins, lipid molecules, water molecules, and ions is constructed using the CHARMM-GUI online platform, and then made electrically neutral using an ion addition algorithm.

[0022] S2 uses multi-stage energy minimization and ensemble equilibrium conditions to perform dynamic simulations of the tactile conformation of the system and obtains the equilibrium system. S21 performs energy minimization on the established system to obtain an energy-stable conformation; Specifically, the three-dimensional biomimetic structural model is subjected to energy minimization processing to eliminate unreasonable interatomic contacts and obtain an energy-stable conformation; S22 used the NVT ensemble to perform temperature equilibration on the energy-stable conformation of the complex, and obtained a temperature-stable simulation system. Specifically, the energy-stable conformation is temperature-balanced using the NVT ensemble to stabilize the system temperature to the target value, thus obtaining a temperature-stable simulation system.

[0023] S23 uses the NVT ensemble to perform pressure balancing on the temperature stability simulation system, thus obtaining the equilibrium simulation system.

[0024] Specifically, the NPT ensemble is used to perform pressure equilibration on the temperature stability simulation system to stabilize the system density and pressure, thus obtaining an equilibrium system for formal simulation.

[0025] S3 combines Fourier transform and linear response theory with the obtained equilibrium system to obtain the vibration spectrum of the desired binding site; The vibration spectrum calculation mentioned refers to establishing the autocorrelation function of the total dipole moment of the entire system and performing Fourier transform and linear response theory to obtain the vibration spectrum of the system. Combined with VMD, the functional groups corresponding to the spectral peaks can be obtained in real time, and the characteristic vibration frequencies corresponding to the ion channel binding sites can be obtained.

[0026] Specifically, this further includes: Calculation of the total dipole moment autocorrelation function. Data on the time-varying total dipole moment of the entire system (including ion channel proteins, the phospholipid bilayer, and the solution environment) under no external field conditions are extracted from the molecular dynamics simulation trajectory of the equilibrium system. Subsequently, the autocorrelation function (ACF) of this total dipole moment is calculated to characterize the intrinsic collective vibrational modes of the system. Vibrational spectrum generation and characteristic frequency extraction. The autocorrelation function obtained in step S31 is subjected to a Fourier transform, converting it from the time domain to the frequency domain, thereby obtaining the vibrational spectrum of the system in the terahertz band (typically 1-20 THz). This spectrum is analyzed to identify significant characteristic peaks; the frequencies corresponding to these peaks are considered key frequencies that may regulate ion channel function. Functional group assignment of characteristic frequencies. Functional group assignment analysis is performed on the characteristic frequencies identified in step S32 using molecular visualization software (such as VMD). By comparing spectral peaks with known vibrational frequencies of specific chemical groups (such as the -C=O stretching vibration of carboxyl groups) in key functional regions of ion channels (e.g., selective filters, SF), the microscopic vibrational sources corresponding to these characteristic frequencies are identified. For example, studies show that a frequency of 51.87 THz may be related to the vibration of carbonyl groups in potassium ion channel selective filters, while a frequency of 42.55 THz may be related to the stretching mode of carboxylate groups in calcium ion channel selective filters. Linear response theory verification is then performed. Based on linear response theory, the resonance condition between the frequency of the external field (terahertz wave) and the intrinsic vibrational characteristic frequency of the system is confirmed. This theoretical framework helps to clarify that when the applied terahertz field frequency matches the characteristic frequency identified in step S33, it can most effectively influence the conformation or ion transport behavior of the ion channel.

[0027] Based on the characteristic frequency derived from the spectrum, S4 loads the corresponding frequency terahertz wave and performs molecular dynamics simulation again. The simulation data is then analyzed in multiple dimensions to obtain quantitative indicators of the interaction between the terahertz wave and the ion channel.

[0028] S41 calculates the root mean square deviation of the conformation trajectory data to obtain the structural stability index of the complex. S42 calculates the change in potential energy before and after loading terahertz waves based on average potential energy. S43 uses an energy decomposition algorithm to analyze the interaction energy of conformation trajectory data and obtains the influence of terahertz waves on the system energy. S44 calculates the effect of terahertz waves on the pore size of ion channels based on aperture calculations. S45 performs dynamic structural analysis on conformational trajectory data based on secondary structure to obtain the influence of terahertz waves on the changes in the structure.

[0029] Specifically, the root mean square deviation of the conformational trajectory data before and after terahertz wave loading was calculated to obtain the system structure stability index; the change of free energy on the key reaction coordinates before and after terahertz wave loading was analyzed based on the mean potential calculation; the influence of terahertz wave on the interaction energy between ion channel and phospholipid and solution was analyzed using the energy decomposition algorithm; the dynamic influence of terahertz wave on the pore size of ion channel was quantified through the pore size analysis algorithm; and the dynamic changes and retention rate of the secondary structure of ion channel under the action of terahertz wave were analyzed based on the secondary structure allocation algorithm.

[0030] Example This invention is based on a molecular dynamics simulation method for the interaction between terahertz waves and potassium ion channels.

[0031] like Figure 1 As shown, the protein-phospholipid system was first constructed. The 7CR0 crystal structure was downloaded from the PDB database, and the portion inserted into the phospholipid bilayer was extracted and saved using VMD software. The prepared protein file was imported, and the membrane system was constructed using CharMM-GUI. Alignment along the z-axis was selected, and the membrane area that could cover the protein was calculated. A POPC membrane was selected and created as needed. Here, a 0.15 mmol / L KCl solution was added as the solution environment, and the membrane protein system was constructed online. Finally, the initial membrane protein structure model's PDB and GRO files were obtained.

[0032] Then, perform dynamic simulations. Prepare the PDB file and topology file, as well as the energy minimization file em.mdp. Generate the energy minimization tpr file using the gmxgrompp command and perform energy minimization. Bond length and bond angle constraints should be maintained at this stage.

[0033] Next, perform restricted dynamics step by step, gradually reducing the restrictions. First, prepare the restricted dynamics .tpr file: run the restricted dynamics simulation, using a smaller 1fs step size for complex systems to avoid simulation crashes at the beginning. Then perform NVT and NPT dynamics simulations.

[0034] Perform a 40ns routine dynamics simulation, which involves preparing a routine dynamics tpr file and running the simulation.

[0035] Finally, a multi-dimensional analysis was conducted.

[0036] First, vibrational spectroscopy calculations are performed to obtain characteristic frequencies. The amino acids to be calculated are reassigned to a group within the ndx file and named, the restrictions on the calculation group are removed, and restrictions are applied to other structures. A spectral calculation file, IR.mdp, is created, with output every four steps, the simulation duration adjusted to 100 ps, ​​and the trajectory and velocity files preserved. Charge information is extracted from the velocity files, and an autocorrelation function is established and a Fourier transform is performed to convert it from the time domain to the frequency domain, thus obtaining the vibrational spectrum of the system in the terahertz band. This calculated spectrum is then combined with VMD software. VMD software can be used to determine the vibrational characteristics of specific chemical bonds, and the required vibrational frequencies at key sites can be obtained through comparison. Figure 2 The frequency here is 51.2 THz.

[0037] The mean potential (PMF) is calculated. An irrelevant ion is selected (in this case, chloride ion as the reference point), and then the interacting ion passing through the ion channel is selected. Both ions are frozen using the `freezegrps` command. The positions of the moving ions are adjusted, and energy minimization is performed. A `pull.mdp` file is created as the potassium ion passage file. Considering the ion passage direction as inward-outward, the traction method is set to `direction`. A dynamic simulation is run to obtain the trajectory `pull.xtc` file. Based on the trajectory file, a `gro` file for each frame position is extracted. A script is used to calculate and extract the number of frames corresponding to each 0.2 Å traction. Finally, umbrella sampling is performed on the extracted frames, and dynamic simulations are conducted for both cases with and without electromagnetic fields. Finally, the `gmxwham` command is used to integrate the calculated results of each frame to obtain the PMF image. Figure 3 As shown, the PMF decreases after loading a specific terahertz wave.

[0038] Calculate the binding energy. Separate the protein and channel ions into two groups and save them to ndx files. Simultaneously, create separate top files for the protein and ions. Process the kinetic XTC and GRO files, using the `gmxtrjcov` command to extract three types of XTC and GRO files: one containing both proteins and ions, one containing only proteins, and one containing only ions. Rerun the kinetic simulation using the `rerun` command, and run `gmxenergy` to calculate the energy for each of the three cases. Subtract the energy of the protein and ion cases from the total total energy to obtain the binding energy. Figure 4 Specific terahertz waves can affect the binding energy of proteins and ions.

[0039] Pore ​​size analysis was performed. The protein portion of the kinetic analysis file was extracted and a specific frame was saved as a PDB file. Cinema software was used to correct erroneous hydrogen bonds in the PDB file before importing it into VMD. Pore size was calculated using the molecular.exe program, with the program's built-in cpoint option set to points within the protein channels, resulting in a dot.vmd_plot file. This file was then loaded into VMD using the source command, allowing visualization of the channel radius. Multiple calculations were performed, and the channel_raddii.dat file was analyzed to obtain pore size changes. Figure 5 Certain terahertz waves can cause changes in the pore size of protein channels.

[0040] Secondary structure analysis is performed. Gromacs 2023 and later versions can directly obtain the changes in each secondary structure using the `dssp` command. Here, we use Gromacs 2018.8 with the `do_dssp` command embedded for secondary structure calculation, employing `gmxdo_dssp`. The command can generate a secondary structure xpm distribution diagram and quantity statistics based on the xtc file. After calculating and statistically analyzing the results, a secondary structure chart can be obtained, such as... Figure 6 Specific terahertz waves can cause changes in structures such as A-helix.

[0041] The above-disclosed embodiments are merely preferred embodiments of a molecular dynamics simulation method based on the interaction between terahertz waves and ion channels according to the present invention. Of course, they should not be construed as limiting the scope of the present invention. Those skilled in the art will understand that all or part of the processes of the above embodiments can be implemented, and equivalent changes made in accordance with the claims of the present invention are still within the scope of the invention.

Claims

1. A molecular dynamics simulation method based on the interaction between terahertz waves and ion channels, characterized in that, Includes the following steps: The target system was modeled to obtain a three-dimensional biomimetic structure of ion channel-phospholipid bilayer-solution environment; Based on the multi-stage energy minimization and ensemble equilibrium conditions, dynamic simulations of the tactile conformation of the system are performed to obtain the equilibrium system. By combining Fourier transform and linear response theory with the obtained equilibrium system, the vibration spectrum of the desired binding site is obtained; Based on the characteristic frequencies derived from the spectrum, terahertz waves of the corresponding frequencies are loaded to perform molecular dynamics simulations again. The simulation data are then analyzed in multiple dimensions to obtain quantitative indicators of the interaction between terahertz waves and ion channels.

2. The molecular dynamics simulation method based on the interaction between terahertz waves and ion channels as described in claim 1, characterized in that, In "Modeling the target system to obtain the three-dimensional biomimetic structure of ion channel-phospholipid bilayer-solution environment", the electroneutrality model can be established by CHARMM-GUI according to the required ion channel type, phospholipid type and liquid environment.

3. The molecular dynamics simulation method based on the interaction between terahertz waves and ion channels as described in claim 1, characterized in that, The process of "dynamic simulation of the thixotropic conformation of the system based on multi-stage energy minimization and ensemble equilibrium conditions to obtain the equilibrium system" includes the following steps: By performing energy minimization on the established system, an energy-stable conformation is obtained; Using the NVT ensemble, the energy-stable conformation of the complex was equilibrated at temperature to obtain a temperature-stable simulation system; Using the NVT ensemble, pressure equilibrium was achieved in the temperature stability simulation system to obtain the equilibrium simulation system.

4. The molecular dynamics simulation method based on the interaction between terahertz waves and ion channels as described in claim 1, characterized in that, In the section "By combining Fourier transform and linear response theory with the obtained equilibrium system, the vibration spectrum calculation mentioned refers to establishing the autocorrelation function of the total dipole moment of the entire system and performing Fourier transform and linear response theory to obtain the vibration spectrum of the system. By combining VMD, the functional groups corresponding to the spectral peaks can be obtained in real time, and the characteristic vibration frequencies corresponding to the ion channel binding sites can be obtained.

5. The molecular dynamics simulation method based on the interaction between terahertz waves and ion channels as described in claim 1, characterized in that, In the section "Based on the characteristic frequencies derived from the spectrum, a corresponding terahertz wave is applied to perform molecular dynamics simulations again, and the simulation data is analyzed in a multi-dimensional molecular manner to obtain quantitative indicators of the interaction between the terahertz wave and the ion channel," the quantitative indicators include: The root mean square deviation of the conformational trajectory data is calculated to obtain the structural stability index of the complex. Based on the average potential calculation, the change in potential energy before and after the application of terahertz waves is obtained; The interaction energy of conformation trajectory data was analyzed using an energy decomposition algorithm to obtain the influence of terahertz waves on the system energy. Based on aperture calculations, the influence of terahertz waves on the aperture size of ion channels was obtained; Dynamic structural analysis of conformational trajectory data based on secondary structure reveals the influence of terahertz waves on the structure.

Images