A molecular dynamics-based method for simulating the rhodopsin photocycle

Abstract

The present invention claims a molecular dynamics rhodopsin protein light cycle simulation research method, which uses monomers for simulation research and belongs to the fields of neuroscience and molecular dynamics. The specific steps are as follows: (1) Use Gaussian software to twist the chromophore C12-C13=C14-C15 dihedral angle and adjust the angle; (2) Embed the monomer protein in the cell membrane environment in CHARMM-GUI software, and set a separate force field for the chromophore small molecule (3) Use NAMD software to perform equilibrium molecular dynamics simulation (4) Analyze the molecular dynamics trajectory, and perform corresponding analysis through the cpptraj module of AMBER and the tcl file of VMD. (5) Calculate free energy using stretching molecular dynamics and umbrella sampling. The inventive method uses molecular dynamics to study the structural changes of rhodopsin during the light cycle and further realizes the simulation of the light cycle, which plays an important role in the clinical application of rhodopsin in the field of neuroscience, brain function research and engineering material design.

Description

Technical Field

[0001] The invention belongs to the fields of neuroscience and molecular dynamics, and in particular relates to a rhodopsin light cycle simulation system based on molecular dynamics. Background Art

[0002] The use of light-sensitive proteins to observe and control cell activities is considered one of the most groundbreaking innovations in the field of neuroscience in recent years. This method has gradually replaced traditional methods for studying neural function. In 2006, Karl Deisseroth first proposed the term optogenetics. Optogenetics is a combination of genetics and optics that can stimulate or inhibit specific events in specific cells of living tissues or animals with behavioral abilities. When isolating the neural circuits of specific neuronal populations, the reversible silencing of neurons plays an important role in their dynamics and behavior. Light-induced functions are controlled by photoreceptor proteins, which contain an embedded chromophore molecule that can sense light. The protein environment promotes the rapid and efficient light response of the chromophore, and the altered chromophore-protein interaction produces changes in the protein structure. The details of protein structural changes must be analyzed at the atomic level. Since the light cycle process can be completed in the range of picoseconds to seconds, the crystal structure of the light cycle intermediates cannot be obtained. Due to the lack of fully open state and intermediate structure crystals of wild-type light-gated ion channels and their mutants, analyzing the structural changes of protein atomic details is a challenge in biophysics and physical chemistry.

[0003] At present, the main research methods for studying the structural information of photosensitive proteins include nuclear magnetic resonance technology, electron microscopy technology, vibrational spectroscopy (such as Raman spectroscopy and infrared spectroscopy) technology, Fourier transform infrared spectroscopy (FTIR) and molecular dynamics technology. The crystal structure obtained by X-ray diffraction can only characterize the instantaneous conformation of the protein, which is often not the active state conformation. Molecular dynamics simulation can construct a protein solvation model that is close to the physiological state, and can study the relationship between the dynamic experimental data established by the computer and the static molecular conformation, and provide dynamic structural image change data that cannot be detected by experiments. Molecular dynamics simulation (MD) is aimed at a multi-body system composed of atomic nuclei and electrons. The motion process of atomic nuclei is simulated by computer, and the conformation and properties of the system are obtained thereby. Any atomic nucleus is regarded as moving according to Newton's law under the potential field provided by other atomic nuclei and electrons in the system. It connects the properties of the microscopic system with the macroscopic system in a computer-assisted way. The photocycle reaction triggered by light absorption leads to proton transport on the visual protein membrane. The C12-C13=C14-C15 dihedral angle and protonation state of the chromophore small molecule in rhodopsin are particularly important for the structure of the photocycle intermediates of opsin. Therefore, improving the chromophore small molecule simulation system plays an important role in studying the molecular mechanism of opsin.

[0004] In the current method of using molecular dynamics to study the structure of opsin, the different protonation states of different intermediates of chromophore small molecules and the C12-C13=C14-C15 dihedral angle model need to be improved. Classical molecular dynamics cannot accurately obtain the conformation of ion channels. There is a huge difference between the time scale of conventional molecular dynamics simulation and the actual biological process, resulting in insufficient sampling. Calculating free energy by stretching molecular dynamics and umbrella sampling can more accurately calculate the free energy change and more accurately analyze the structural mechanism of ion channels. By improving the force field environment of chromophore small molecules and combining the method of stretching molecular dynamics to calculate free energy, the structural characteristics of gated channel proteins can be more accurately obtained. Understanding the structural characteristics of opsin and further improving the structural characteristics of opsin are conducive to application in basic research applications and clinical applications, creating possibilities for similar gene therapy for vision loss or loss caused by retinal diseases, cardiac conduction system diseases and nervous system diseases. Summary of the invention

[0005] The present invention aims to solve the above problems of the prior art. A method for simulating the photocycle of rhodopsin based on molecular dynamics is proposed. The technical solution of the present invention is as follows:

[0006] A method for simulating the rhodopsin photocycle based on molecular dynamics comprises the following steps:

[0007] Step 1: Use Gaussian software to twist and adjust the dihedral angle of C12-C13=C14-C15 of the monomer PDB file chromophore small molecule;

[0008] Step 2: embed the monomeric protein after the dihedral angle of the chromophore small molecule is twisted into the cell membrane environment in the CHARMM-GUI software, set a separate force field for the chromophore small molecule, and prepare the molecular dynamics simulation system;

[0009] Step 3: Use NAMD software to perform equilibrium molecular dynamics simulation, extract representative conformations of protein structure, further twist the C12-C13=C14-C15 dihedral angle of small molecule, fix the C12-C13=C14-C15 dihedral angle, perform molecular dynamics simulation again until the dihedral angle bond angle twist is completed, and perform molecular dynamics simulation without fixing the dihedral angle of small molecule;

[0010] Step 4: Analyze the trajectory after molecular dynamics simulation through the cpptraj module of AMBER and the tcl file of VMD;

[0011] Step 5: Extract representative conformations from molecular dynamics trajectories and calculate free energy using stretched molecular dynamics and umbrella sampling.

[0012] Furthermore, the step 1, using Gaussian software to twist the chromophore C12-C13=C14-C15 dihedral angle and adjust the angle, specifically includes: manually deleting other polymers from the protein multimer file downloaded from the PDB database to leave only monomer crystals, retaining single-chain amino acid residues, deleting all other crystal water molecules, and loading the monomer PDB file using Gaussian software; selecting the C12-C13=C14-C15 dihedral angle double bond of the isomerized chromophore small molecule to isomerize by rotating it by 20°, and keeping it fixed in each torsion equilibrium simulation until the C12-C13=C14-C15 dihedral angle reaches about 0°.

[0013] Furthermore, in step 2, the monomer protein is embedded in the cell membrane environment in the CHARMM-GUI software, and missing hydrogen atoms and residues are added; the monomer is embedded in a periodic box with POPC, TIP, and 0.15 mol·L-1 NaCl membrane on the CHARMM-GUI platform, and the entire monomer is loaded into the force field of CHARMM36m;

[0014] Furthermore, in the step 2, for the simulation system where the chromophore needs to be isomerized, during the energy optimization process, the monomer is kept fixed and the C12-C13=C14-C15 dihedral angle of the chromophore is fixed, and the system is energy minimized; when the dihedral angle reaches the isomerization angle to complete the required angle of 0°, and the C12-C13=C14-C15 dihedral angle is not fixed, molecular dynamics energy minimization and equilibrium dynamic simulation are performed;

[0015] The chromophore is also set as a deprotonated simulation system, and the deprotonated chromophore is placed in a separate force field. This ensures that the chromophore small molecule will not be broken in the force field of CHARMM36m.

[0016] Furthermore, the molecular dynamics energy minimization and equilibrium dynamics simulation are performed for the case where the C12-C13=C14-C15 dihedral angle is not fixed; during the energy minimization process, the protein is fixed and only water molecules are allowed to move, so that the solvent in the system is optimized. The equilibrium dynamics simulation is simulated under the canonical ensemble (NVT), by adjusting the speed of the atoms to keep the system kinetic energy constant, relaxing the temperature T to a steady state, and then switching to constant-pressure, constant-temperature (NPT) for pre-equilibrium, and by adjusting the volume of the system to keep the pressure constant, so that the pressure P is relaxed to a steady state. After that, the system can perform actual simulations.

[0017] Furthermore, the step three, using NAMD software to perform equilibrium molecular dynamics simulation, specifically includes: the first stage of the equilibrium step is carried out under the canonical ensemble (NVT), that is, the system temperature, volume and particle number remain unchanged, the system is slowly heated under the NVT ensemble, and at the required temperature, the initial velocities of all atoms are randomly distributed in the Maxwell-Boltzmann distribution; finally, the protein system is allowed to run the dynamics for a long enough time under isothermal and isobaric (constant-pressure, constant-temperature, NPT) with constant temperature, pressure and particle number so that the protein structure tends to be stable. This process will adjust the pressure of the system and thus converge the system density.

[0018] Furthermore, the step 4 is to analyze the molecular dynamics trajectory, and the corresponding analysis is performed through the cpptraj module of AMBER and the tcl file of VMD, specifically including: the atoms in the system move within a time step according to the initial speed, and then the new coordinates of the atoms are recorded in the trajectory file; after a time step, the atoms continue to move, and then the potential energy and speed of the atoms are recalculated and recorded in the trajectory file; and this process is repeated until the entire simulation is completed; after the protein motion trajectory tends to be stable, the stable motion trajectory is used to extract data for analyzing the molecular mechanism of the protein; the cpptraj functional module of AMBER is used to analyze the representative structural characteristics of rhodopsin and the distribution of water molecules and hydrogen bonds in the ion channel; the VMD software is used to calculate the number of water molecules in the ion channel to indirectly reflect the degree of helical hydration of the ion channel.

[0019] Furthermore, the step 5, step 5, uses stretched molecular dynamics and umbrella sampling to calculate free energy, specifically including:

[0020] The cluster representative conformations obtained in the equilibrium state are used as models, and the free energy is calculated by using an umbrella sampling method after the ions are stretched from one side of the ion channel to the other side at a constant speed. In the SMD research method, the pulled atoms are connected to an imaginary spring and an imaginary atom, and the magnitude of the pulling force is determined by the deformation and spring coefficient of the spring, which satisfies Hooke's law:

[0021] F(t)=2k(vt-x(t))

[0022] k is the elastic constant of the spring, v is the velocity of the pulled atom, and x(t) is the position of the pulled atom at time t.

[0023] Umbrella sampling uses WHAM, Weighted Histogram Analysis Method, with the aim of converting the statistical results of biased sampling into those of unbiased sampling. The basic idea is to first generate several conformations of the system along the reaction coordinate (z) at a certain interval, and each conformation corresponds to a window i. A bias potential Wi(Z) (generally a simple harmonic spring potential) is added to each umbrella sampling window to overcome the problem of insufficient sampling of certain energetically unfavorable conformations. Then, an independent bias simulation is performed on each umbrella sampling window, and the umbrella histograms in each sampling window are combined to represent the probability distribution of the conformation along the reaction coordinate under the action of the bias potential.

[0024] The basic principle of calculating the mean force potential by umbrella sampling is:

[0025]

[0026] The advantages and beneficial effects of the present invention are as follows:

[0027] (1) The present invention sets up a more reasonable simulation system for the isomerization model of the chromophore small molecule that plays an important role in the photocycle of rhodopsin. Compared with the previous simulation system, a separate chromophore force field and a dihedral angle fixed method are set up to perform the energy minimization molecular dynamics model.

[0028] (2) During the equilibrium trajectory analysis, the water molecule analysis and hydrogen bond analysis, which are important for the openness of rhodopsin, were completed using the cpptraj module of AMBER. In addition, the statistical method for the number of water molecules in the rhodopsin ion channel was implemented using the tcl script of VMD.

[0029] (3) A model system for free energy calculation was established. The free energy was calculated by combining the stretching molecular dynamics with constant stretching and the umbrella sampling method. This fully improved the shortcomings of insufficient sampling in classical molecular dynamics and helped to understand the ionic conductivity of rhodopsin by understanding the energy that needs to be overcome during the permeation process. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 The present invention provides a preferred embodiment Figure 1 This is the dihedral angle change model of the chromophore small molecule C12-C13=C14-C15 of the present invention;

[0031] Figure 2 It is the key code of the force field of the chromophore small molecule deprotonation model of the present invention;

[0032] Figure 3 The protein monomer periodic box model of the present invention;

[0033] Figure 4 The key code for energy minimization of the fixed chromophore small molecule C12-C13=C14-C15 of the present invention;

[0034] Figure 5 It is the cluster analysis model of the present invention;

[0035] Figure 6 is the hydrogen bond analysis model of the present invention;

[0036] Figure 7 The water distribution analysis model of the present invention;

[0037] Figure 8 The water molecule quantity analysis model of the present invention;

[0038] Fig. 9 The free energy calculation and analysis model of the present invention;

[0039] Fig.10The figure is a flow chart of the rhodopsin photocycle simulation method based on molecular dynamics of the present invention. DETAILED DESCRIPTION

[0040] The following will describe the technical solutions in the embodiments of the present invention in detail in conjunction with the accompanying drawings in the embodiments of the present invention. The described embodiments are only part of the embodiments of the present invention.

[0041] The technical solution of the present invention to solve the above technical problems is:

[0042] To achieve the above object, the present invention adopts the following technical solution: a rhodopsin light cycle simulation system based on molecular dynamics, comprising the following steps:

[0043] Step 1: Use Gaussian software to twist the chromophore C12-C13=C14-C15 dihedral angle and adjust the angle; manually delete the redundant protein multimer file downloaded from the PDB database to turn it into a monomer crystal, retain the single-chain amino acid residues, delete all other small molecules and other structures. Use Gaussian software to load the monomer PDB file. Select the structure of the chromophore small molecule isomerization C12-C13=C14-C15 dihedral angle double bond to isomerize with a 20° rotation and keep it fixed in each torsion equilibrium simulation until the C12-C13=C14-C15 dihedral angle reaches approximately 0°.

[0044] Step 2: Embed the monomer protein into the cell membrane environment in the CHARMM-GUI software, add missing hydrogen atoms and residues; embed the monomer into a periodic box with POPC, TIP, and 0.15 mol·L-1 NaCl membrane on the CHARMM-GUI platform, and load the entire monomer into the force field of CHARMM36m.

[0045] 1. For the simulation system where the chromophore needs to be isomerized, during the energy optimization process, the monomer and the C12-C13=C14-C15 dihedral angle of the chromophore are kept fixed, and the energy of the system is minimized. When the dihedral angle reaches the required angle of 0° for isomerization, the molecular dynamics energy minimization and equilibrium dynamic simulation are performed with the C12-C13=C14-C15 dihedral angle not fixed.

[0046] 2. The chromophore needs to be set as a deprotonated simulation system, and the deprotonated chromophore should be placed in a separate force field. This can ensure that the chromophore small molecule will not be broken in the force field of CHARMM36m.

[0047] Step 3: Equilibrium molecular dynamics simulation; The first stage of the equilibrium step is carried out under the NVT ensemble, that is, the system temperature, volume and number of particles remain unchanged. The system is slowly heated up under the NVT ensemble. Usually at the required temperature, the initial velocities of all atoms can be randomly distributed in the Maxwell-Boltzmann distribution. Finally, the protein system is allowed to run the dynamics in the NPT ensemble with constant temperature, pressure and number of particles for a long enough time to stabilize the protein structure. This process will adjust the pressure of the system and converge the system density. After completing these equilibrium steps, the system can be actually simulated.

[0048] Step 4: Analyze the molecular dynamics trajectory through the AMBER cpptraj module and VMD tcl file. The atoms in the system move within a time step according to the initial velocity, and then the new coordinates of the atoms are recorded in the trajectory file; after a time step, the atoms continue to move, and then the potential energy and velocity of the atoms are recalculated and recorded in the trajectory file; and this process is repeated until the entire simulation is completed. After the protein motion trajectory tends to be stable, the stable motion trajectory is used to extract data for analyzing the molecular mechanism of the protein. The AMBER cpptraj function module can be used to analyze the representative structural features of rhodopsin and the distribution of water molecules and hydrogen bonds in ion channels. The VMD software can be used to calculate the number of water molecules in the ion channel to indirectly reflect the degree of helical hydration of the ion channel.

[0049] Step 5. Calculate the free energy simulation. Classical molecular dynamics easily causes sampling to fall into the energy minimum in the conformational space, making it difficult to realize the process of ion transmembrane conduction. In order to improve the problem of insufficient spatial sampling caused by the huge difference in time scale between classical molecular dynamics simulation and real biological processes. The method of combining stretching molecular dynamics and umbrella sampling with calculating free energy can solve this problem well. In order to better study the thermodynamic properties of ion transmembrane conduction, the constant rate stretching method can better study the changes in the internal structure of proteins. The umbrella sampling simulation can find the equilibrium state in a set of umbrella sampling windows, and obtain a relatively accurate PMF curve by re-weighting the constant bias potential well along the reaction coordinate. This provides a good theoretical basis for subsequent optogenetic gene modification.

[0050] First, a dimer crystal named 6EDQ was downloaded from the PDB database, leaving only the monomer structure, and Gaussian software was used to twist the C12-C13=C14-C15 dihedral angle of the small molecule by 20°.

[0051] Figure 1 、Structure model of chromophore small molecule. The initial angle is -179.84252, and the C12-C13=C14-C15 angle is rotated 20° each time.

[0052] The rotated monomer was structurally optimized in Membrane Builder, including adding side chains of amino acid residues, hydrogen atoms, assigning protonation states, and optimizing amino acid side chains. Finally, we embedded the optimized structure into the phospholipid bilayer (POPC), placed the system in the cubic periodic box of the TIP3P solvent model, and added counterions (0.15 mol·L-1 NaCl) to maintain the electrical neutrality of the system. The entire monomer was loaded into the force field of CHARMM36m. In this step, the chromophore small molecule Schiff base can be deprotonated and the protonation states of other acidic and basic amino acids can be checked. The separate force field prm file and rtf file for the deprotonated Schiff base can be loaded. See the key code for details. Figure 2 :

[0053] Figure 2 、Key codes of small molecule deprotonation force field. A. Key codes of prm file (topology file of deprotonated retinal). B. Key codes of rtf file (retinal related parameters).

[0054] After the chromophore small molecule structure and the required force field file are loaded, its periodic box model is as follows Figure 3 :After obtaining the initial model system, in the simulation parameter setting, the PME (Particle Mesh Ewald) algorithm is used to deal with long-range electrostatic interactions, and the cutoff value of short-range non-bonded interactions is set to The SHAKE algorithm is used to restrict the vibration of the covalent bonds connected to all hydrogen atoms, and the step size is set to 2fs. The C12-C13=C14-C15 dihedral angle is fixed, and the system is minimized by the steepest descent method. Then the system is heated in the NVT ensemble from 0K to 310K. During the heating process, the Langevin piston algorithm is used to keep the pressure at 1 bar. The key code for fixing the C12-C13=C14-C15 dihedral angle is as follows Figure 4 :

[0055] Then, the NPT ensemble was used to perform molecular dynamics simulations without any restrictions on the system under constant temperature and pressure conditions. Langevin dynamics was used during the simulation to ensure that the temperature in all simulations was constant at 310K, and the damping coefficient was taken as 1ps-1. The system pressure was maintained at a constant pressure of 1 atmosphere using the Nose-Hoover Langevin piston method, with a decay period of 25fs and a damping period of 50fs. A relatively stable trajectory file of the protein was obtained, and cluster analysis was performed using the cpptraj module of the AMBER software. The clustering comparison of the chromophore small molecule of GtACR1 before and after isomerization is shown in the figure. Figure 5 :

[0056] The degree of ion development of rhodopsin is closely related to the hydrogen bond network and the degree of helical hydration. Therefore, it is important to understand the hydrogen bond network and water molecule distribution of rhodopsin. Figure 6 For water density analysis using cpptraj, Figure 7 This is a hydrogen bond analysis. The number of hydrogen bonds in the ion channel and the distribution of water molecules can provide insights into the openness of the ion channel.

[0057] Although the distribution of water molecules has a certain influence on the openness of ion channels, being able to calculate the number of water molecules in ion channels can also indirectly reflect the openness and hydration degree of ion channels. The number of water molecules in the range, such as Figure 8 :

[0058] The representative conformation of classical molecular dynamics is used as the initial conformation of stretched molecular dynamics. In the simulation, an imaginary external force is artificially applied to Cl-, so that the ion can move along the channel. Then, detailed information about the atomic coordinates during the entire stretching process is obtained, and important information about the process of ion transmembrane conduction is obtained by analyzing the changes in structure and energy. During the stretching process, a hard spring is used to pull the SMD atoms at a constant speed. In order to avoid the translation and rotation of the protein in the SMD simulation, an amino acid residue outside the seven helices of the protein is added. The reaction coordinates used in the umbrella sampling simulations are the same as those used in the SMD simulations, i.e., the z distance along the channel axis. In order to obtain fine structural information, the reaction coordinate space is Divide a window. Use the SMD trajectory to extract the structure near the center of each window as the initial conformation of each independent window, and then constrain the ion in the window center in the z direction. Finally, the weighted histogram analysis method (WHAM) is used to construct the PMF based on these umbrella sampling simulation trajectories, such as Fig. 9 .

[0059] The systems, devices, modules or units described in the above embodiments may be implemented by computer chips or entities, or by products with certain functions. A typical implementation device is a computer. Specifically, the computer may be, for example, a personal computer, a laptop computer, a cellular phone, a camera phone, a smart phone, a personal digital assistant, a media player, a navigation device, an email device, a game console, a tablet computer, a wearable device, or a combination of any of these devices.

[0060] It should also be noted that the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, commodity or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, commodity or device. In the absence of more restrictions, the elements defined by the sentence "comprises a ..." do not exclude the existence of other identical elements in the process, method, commodity or device including the elements.

[0061] The above embodiments should be understood to be only used to illustrate the present invention and not to limit the protection scope of the present invention. After reading the contents of the present invention, technicians can make various changes or modifications to the present invention, and these equivalent changes and modifications also fall within the scope defined by the claims of the present invention.

Claims

1. A method for simulating the photocycle of rhodopsin based on molecular dynamics, characterized in that: The following steps are involved: Step 1: Use Gaussian software to twist and adjust the dihedral angle of C12-C13=C14-C15 of the monomer PDB file chromophore small molecule; Step 2: embed the monomeric protein after the dihedral angle of the chromophore small molecule is twisted into the cell membrane environment in the CHARMM-GUI software, set a separate force field for the chromophore small molecule, and prepare the molecular dynamics simulation system; Step 3: Use NAMD software to perform equilibrium molecular dynamics simulation, extract representative conformations of protein structure, further twist the C12-C13=C14-C15 dihedral angle of small molecule, fix the C12-C13=C14-C15 dihedral angle, perform molecular dynamics simulation again until the dihedral angle bond angle twist is completed, and perform molecular dynamics simulation without fixing the dihedral angle of small molecule; Step 4: Analyze the trajectory after molecular dynamics simulation through the cpptraj module of AMBER and the tcl file of VMD; Step 5: Extract representative conformations from molecular dynamics trajectories and calculate free energy using stretched molecular dynamics and umbrella sampling. The second step is to embed the monomer protein into the cell membrane environment in the CHARMM-GUI software, add missing hydrogen atoms and residues; embed the monomer into a periodic box of the monomer coated with POPC, TIP, and 0.15 mol·L-1 NaCl membrane on the CHARMM-GUI platform, and load the entire monomer into the force field of CHARMM36m; In the step 2, for the simulation system where the chromophore needs to be isomerized, during the energy optimization process, the monomer is kept fixed and the C12-C13=C14-C15 dihedral angle of the chromophore is fixed, and the system is energy minimized; when the dihedral angle reaches the isomerization angle to complete the required angle of 0°, the molecular dynamics energy minimization and equilibrium state dynamic simulation are performed for the case where the C12-C13=C14-C15 dihedral angle is not fixed; For the simulation system where the chromophore is set to be deprotonated, the deprotonated chromophore is placed in a separate force field. This ensures that the chromophore small molecule will not be broken in the force field of CHARMM36m.

2. The method for simulating the rhodopsin light cycle based on molecular dynamics according to claim 1, characterized in that: The step 1, using Gaussian software to twist the chromophore C12-C13=C14-C15 dihedral angle and adjust the angle, specifically includes: manually deleting other polymers from the protein polymer file downloaded from the PDB database to leave only monomer crystals, retaining single-chain amino acid residues, deleting all other crystal water molecules, and loading the monomer PDB file using Gaussian software; selecting the C12-C13=C14-C15 dihedral angle double bond of the isomerized chromophore small molecule to isomerize by rotating it by 20°, and keeping it fixed in each torsion equilibrium simulation until the C12-C13=C14-C15 dihedral angle reaches approximately 0°.

3. The method for simulating the rhodopsin light cycle based on molecular dynamics according to claim 1, characterized in that: The molecular dynamics energy minimization and equilibrium dynamics simulation are then performed for the case where the C12-C13=C14-C15 dihedral angle is not fixed. During the energy minimization process, the protein is fixed and only water molecules are allowed to move, so that the solvent in the system is optimized; The equilibrium dynamic simulation is simulated under the canonical ensemble NVT. The kinetic energy of the system is kept constant by adjusting the speed of atoms, the temperature T is relaxed to a steady state, and then the isothermal and isobaric NPT is switched for pre-equilibrium. The pressure is kept constant by adjusting the volume of the system, thereby relaxing the pressure P to a steady state.

4. The method for simulating the rhodopsin light cycle based on molecular dynamics according to claim 3, characterized in that: The step three, using NAMD software to perform equilibrium molecular dynamics simulation, specifically includes: the first stage of the equilibrium step is carried out under the canonical ensemble NVT, that is, the system temperature, volume and particle number remain unchanged, the system is slowly heated under the NVT ensemble, and at the required temperature, the initial velocities of all atoms are randomly distributed in the Maxwell-Boltzmann distribution; finally, the protein system is allowed to run the dynamics for a long enough time under isothermal and isobaric NPT with constant temperature, pressure and particle number so that the protein structure tends to be stable. This process will adjust the pressure of the system and thus converge the system density.

5. The method for simulating the rhodopsin light cycle based on molecular dynamics according to claim 3, characterized in that: The fourth step is to analyze the molecular dynamics trajectory, and perform corresponding analysis through the cpptraj module of AMBER and the tcl file of VMD, which specifically includes: the atoms in the system move within a time step according to the initial speed, and then the new coordinates of the atoms are recorded in the trajectory file; after a time step, the atoms continue to move, and then the potential energy and speed of the atoms are recalculated and recorded in the trajectory file; and this process is repeated until the entire simulation is completed; after the protein motion trajectory tends to be stable, the stable motion trajectory is used to extract data for analyzing the molecular mechanism of the protein; the cpptraj function module of AMBER is used to analyze the representative structural characteristics of rhodopsin and the distribution of water molecules and hydrogen bonds in the ion channel; the VMD software is used to calculate the number of water molecules in the ion channel to indirectly reflect the degree of helical hydration of the ion channel.

6. The method for simulating the rhodopsin light cycle based on molecular dynamics according to claim 5, characterized in that: The step 5, calculating the free energy by using stretching molecular dynamics and umbrella sampling, specifically includes: The cluster representative conformations obtained in the equilibrium state are used as models, and the free energy is calculated by using an umbrella sampling method after the ions are stretched from one side of the ion channel to the other side at a constant speed. In the SMD research method, the pulled atoms are connected to an imaginary spring and an imaginary atom, and the magnitude of the pulling force is determined by the deformation and spring coefficient of the spring, which satisfies Hooke's law: F(t)=2k(vt-x(t)) k is the elastic coefficient of the spring, v is the velocity of the pulled atom, and x(t) is the position of the pulled atom at time t; Umbrella sampling uses the WHAM weighted histogram analysis method, the purpose of which is to convert the statistical results of biased sampling into the statistical results of unbiased sampling. The basic idea is to first generate several conformations of the system along the reaction coordinate (z) at a certain interval size, each conformation corresponds to a window i, and add a bias potential W in each umbrella sampling window. i (Z), to overcome the problem of undersampling of some energetically unfavorable conformations, an independent bias simulation is then performed for each umbrella sampling window, and the umbrella histograms in each sampling window are combined to represent the probability distribution of conformations along the reaction coordinate under the bias potential; The basic principle of calculating the mean force potential by umbrella sampling is:

Images