A method for path recognition of ligand-gated ion channels based on molecular dynamics simulation

Abstract

The present invention discloses a method for identifying the path of ligand-gated ion channels based on molecular dynamics simulation. The path identification method includes releasing the ligand-gated ion channel from the constrained state by using conventional molecular dynamics simulation, simulating the random process of ligand binding to the ligand-gated ion channel, and then performing targeted molecular dynamics simulation. After the targeted molecular dynamics simulation, conventional molecular dynamics simulation is performed again to further optimize the structure of the ligand-gated ion channel so that it is in a relaxed state. Gaussian accelerated molecular dynamics simulation is used to smooth the potential energy surface of the ligand-gated ion channel and reduce the energy barrier. Further combined with methods such as dynamic network analysis, the binding sites and allosteric paths of the ligand-gated ion channel and the ligand are analyzed, which is of great significance for studying the mechanism of action of ligand-gated ion channels and ligands binding and developing new allosteric regulators.

Description

Technical Field

[0001] This application relates to the technical fields of computer technology, biomacromolecular multiscale simulation technology, and computational structural biology technology, and particularly relates to a method and device for identifying paths of ligand-gated ion channels based on molecular dynamics simulation. Background Art

[0002] Ligand-gated ion channels (LGICs) are a class of protein complexes widely distributed on cell membranes. Driven by an electrochemical gradient, they mediate the flow of charged ions from the high-concentration side to the low-concentration side. When LGICs bind to specific ligands, they will cause specific conformational changes, thereby forming the open and closed states of the channels. LGICs regulate the opening and closing states of the channels through the allosteric effect generated by ligand binding. Its allosteric communication is a biological process promoted by long-range amino acid communication in space. When a ligand binds or an amino acid at a distant site changes, it will remotely affect the active site, and this effect leads to the opening or closing of the channel through multiple conduction paths.

[0003] Molecular dynamics simulation is a method based on Newtonian mechanics, quantum mechanics, and statistical mechanics, which uses a computer to solve the motion equations of a molecular system. It has been widely used in the study of the allosteric effects of proteins. By analyzing conformational motions and motion correlations, etc., it reveals the allosteric paths and key residues of the protein system.

[0004] The allosteric effect of LGICs usually occurs on a time scale of milliseconds (ms, 10 -3 -3 s), and traditional conventional molecular dynamics simulation (cMD) can generally only capture conformational changes on a time scale of nanoseconds (ns, 10 -9 -9 s) or microseconds (μs, 10 -6 -6 s). Therefore, it may be difficult to capture hidden conformational changes and allosteric paths. There is usually no allosteric pocket in the apo state structure without bound ligands, and only when ligands are present and the allosteric sites are in a relaxed state do they dominate in the conformational ensemble. Therefore, in order to capture hidden conformational changes, it is necessary to first bind ligands to LGICs. After LGICs bind to ligands, they will induce conformational changes in the surrounding domains, and these conformational changes are transmitted to the pore domain through multiple allosteric paths, causing the channel to open.

[0005] In the prior art, the binding pathway between dopamine receptors and ligands was studied by combining kinetic simulations. Twenty-nine molecular dynamics simulations were repeated to screen out the residue sites that play a key role in the binding pathway (Thomas T, Fang Y, Yuriev E, et al. Ligand Binding Pathways of Clozapine and Haloperidol in the Dopamine D2 and D3 Receptors [J]. Journal of Chemical Information & Modeling, 2016: acs.jcim.5b00457.). However, the structure of LGICs is complex, and it is difficult to capture complex allosteric pathways only by repeated kinetic simulations. Due to the large size of the LGICs system and the large time scale of conformational changes, traditional molecular dynamics simulations alone require a large amount of manpower and material resources, and their simulation effects are difficult to cover all conformational changes. Regarding the FMRFamide-activated sodium channel (FaNaC), in the prior art, site-directed mutagenesis, electrophysiology, and molecular dynamics simulations were combined to study its structure and mechanism (Fenglian Liu, Yu Dang, Lu Li, et al. Structure and mechanism of a neuropeptide-activated channel in the ENaC / DEG superfamily, Nature Chemical Biology, 2023, 19(10): 1276-1285), and the residues that play a key role in conformational changes were screened out, but the allosteric binding pathway of the ligand and the receptor was not obtained only by relying on molecular dynamics simulations.

[0006] Molecular docking technology is often used to simulate the binding of ligands to LGICs. Common molecular docking technologies include rigid docking. In this method, during the calculation process, the conformations of the molecules participating in the docking do not change, only the spatial positions and postures of the molecules are changed. The method is simple and the computational amount is relatively small, but it is easy to ignore the flexibility of the molecules, resulting in a decrease in the accuracy of the docking results. Moreover, rigid docking cannot identify the conformational changes that may occur during the binding process of molecules, missing some important binding modes and unable to truly reflect the actual interaction strength. Therefore, random binding, as an important prerequisite, is of great significance for studying the allosteric effects generated by ligand binding to understand the gating mechanism of LGICs and develop new allosteric drugs. Summary of the Invention

[0007] Based on the above, the present invention provides a method for identifying the path of ligand-gated ion channels based on molecular dynamics simulation or a method for molecular dynamics simulation of the binding path and allosteric path between ligand-gated ion channels and ligands. Through various simulation techniques such as conventional molecular dynamics, targeted molecular dynamics, and Gaussian accelerated molecular dynamics, the random binding of ligands to ligand-gated ion channels is simulated to obtain the allosteric path of ligand-gated ion channels from the Apo state to the Holo state. By using analysis methods such as clustering analysis, binding free energy analysis, interaction analysis, and dynamic network analysis, the ligand binding path of ligands entering the binding pocket is analyzed, and then the allosteric path between the ligand and the pore domain is analyzed. The method described in this application realizes the random binding of ligands to ligand-gated ion channels, can obtain the intermediate state structure of ligand-gated channels and the allosteric path between the ligand and the pore domain in a relatively short simulation time, especially identifies more and newer sites than the prior art, and determines the shortest allosteric path, providing a new research idea for studying its gating mechanism and developing new allosteric modulators to regulate this channel.

[0008] In the first aspect of the present application, a method for identifying the path of ligand-gated ion channels based on molecular dynamics simulation is provided. The path identification method includes:

[0009] Step 1: Obtain the crystal structures of the Holo state and Apo state of the ligand-gated ion channel as the initial structures, and construct a simulation system;

[0010] Step 2: Perform conventional molecular dynamics simulation on the simulation system to obtain the initial trajectory of the binding of the ligand-gated ion channel to the ligand and extract the representative pose structure of the binding of the ligand-gated ion channel to the ligand;

[0011] Step 3: Perform targeted molecular dynamics simulation on the representative pose structure of the binding of the ligand-gated ion channel to the ligand and the crystal structure of the Holo state, and extract the initial intermediate state structure;

[0012] Step 4: Perform conventional molecular dynamics simulation on the initial intermediate state structure to obtain the stable trajectory of the intermediate state and extract the stable structure of the intermediate state;

[0013] Step 5: Perform Gaussian accelerated molecular dynamics simulation on the stable structure of the intermediate state to obtain the Gaussian accelerated molecular dynamics simulation trajectory;

[0014] Step 6: Perform dynamic network analysis on the stable trajectory of the intermediate state obtained in Step 4 and the Gaussian accelerated molecular dynamics simulation trajectory obtained in Step 5. Preferably, it also includes correlation coefficient analysis and / or shortest path analysis to obtain the allosteric path of the binding of the ligand to the ligand-gated ion channel.

[0015] The Holo state is the ligand-gated ion channel bound to the ligand, and the Apo state is the ligand-gated ion channel not bound to the ligand.

[0016] In the conventional molecular dynamics simulation in the path recognition method, it includes simulating the random binding of a ligand to a ligand-gated ion channel.

[0017] Specifically, obtaining the crystal structures of the Holo state and the Apo state of the ligand-gated ion channel in step one includes: obtaining the crystal structures of the Holo state and the Apo state of the ligand-gated ion channel from a database, deleting other components in the crystal structure of the Apo state except for the target ligand-gated ion channel, and deleting other components in the crystal structure of the Holo state except for the target ligand-gated ion channel and the ligand; reconstructing the missing structural regions in the crystal structure.

[0018] Preferably, the database includes the PDB database or the RCSB database.

[0019] Preferably, reconstructing the missing structural regions in the crystal structure includes predicting the missing SpecificInsertionⅡ region and then completing it.

[0020] In a specific embodiment of the present invention, reconstructing the missing structural regions in the crystal structure includes predicting the structures of the missing amino acid parts using AlphaFold2 and RoseTTAFold respectively, verifying the optimal structure, then further optimizing the structure using SwissModel homology modeling, and verifying the structural rationality using a Ramachandran plot.

[0021] Specifically, the simulation system in step one includes the crystal structure of the Apo or Holo state of the ligand-gated ion channel; the simulation system also includes a ligand, a phospholipid bilayer, solvent molecules, and ions.

[0022] Preferably, the distance between the ligand-gated ion channel and the ligand is above, such as 15, 16, 17, 18, 19, and above.

[0023] More preferably, the distance between the extracellular domain of the ligand-gated ion channel and the ligand is above, such as 15, 16, 17, 18, 19, and above.

[0024] In a specific embodiment of the present invention, the ligand and the ligand-gated ion channel or its extracellular domain are evenly arranged in a square at a distance not less than to ensure the random binding process of the ligand.

[0025] Preferably, the total number of atoms in the simulation system is 200,000 - 400,000, such as 200,000, 250,000, 300,000, 350,000, 400,000.

[0026] In a specific embodiment of the present invention, step 1 specifically includes:

[0027] Step 1.1: Obtain the crystal structures of the Holo state and Apo state of the ligand-gated ion channel from the PDB database;

[0028] Step 1.2: Delete other components in the crystal structure of the Apo state except the target ligand-gated ion channel, and delete other components in the crystal structure of the Holo state except the target ligand-gated ion channel and the ligand; reconstruct the missing structural regions in the crystal structure;

[0029] Step 1.3: Randomly add ligands to the crystal structure of the Apo state of the ligand-gated ion channel;

[0030] Step 1.4: Construct a simulation system in the CHARMM-GUI software, and the simulation system includes the crystal structure of the Apo state or Holo state of the ligand-gated ion channel; the simulation system also includes ligands, phospholipid bilayers, solvent molecules and ions.

[0031] Preferably, the solvent molecules in step 1.4 include water molecules.

[0032] Step 2 specifically includes:

[0033] Step 2.1: Minimize the energy of the simulation system and heat it, and perform unconstrained conventional molecular dynamics simulation in the NPT ensemble. Preferably, run the simulation system to a relatively equilibrium state;

[0034] Step 2.2: Calculate the root mean square deviation of the ligand heavy atoms, the contact probability between the ligand and the binding site residues, and the binding free energy between the ligand and the ligand-gated ion channel. Preferably, it also includes the interaction between the residues and the ligand, and determine the ligand binding site and / or binding path;

[0035] Step 2.3: Perform clustering analysis on the initial trajectory of the binding of the ligand-gated ion channel to the ligand, and extract the centroid structure of each cluster as the representative pose structure of the binding of the ligand-gated ion channel to the ligand.

[0036] Preferably, the energy minimization of the simulation system described in step 2.1 is achieved by adding restraining forces to the ligand-gated ion channel and the lipid.

[0037] The restraining force is 0.5 - 15 kcal.mol -1 . Any value within the range. Preferably, the restraining force is 1 - 10 kcal.mol -1 . Any value within the range, such as 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 kcal.mol -1 .

[0038] In a specific embodiment of the present invention, by adding a restraining force of 10.0 kcal.mol-1 to the ligand-gated ion channel and 5.0 kcal.mol-1 to the lipid, energy minimization is achieved. of restraining force, adding 5.0 kcal.mol-1 to the lipid. of restraining force to achieve energy minimization.

[0039] Preferably, the heating in step 2.1 includes raising the temperature of the simulation system in the NVT ensemble to 300 - 350 K within 100 - 300 ps. Preferably, the heating includes raising the temperature of the simulation system to 300 - 320 K within 200 - 300 ps.

[0040] In a specific embodiment of the present invention, the heating is to heat the temperature of the simulation system in the NVT ensemble to 310 K within 250 ps.

[0041] The conventional molecular dynamics simulation in step 2.1 is performed in NAMD2 and NAMD3.

[0042] The step size of the conventional molecular dynamics simulation is 1 - 2 fs, and the trajectory is saved every 20 - 100 ps. Preferably, the step size of the conventional molecular dynamics simulation is 2 fs, and the trajectory is saved every 100 ps.

[0043] The conventional molecular dynamics simulation in step 2.1 includes setting the cut-off distance of non-bonded interactions to any value within the range, such as Preferably, the cut-off distance of the non-bonded interactions is

[0044] The conventional molecular dynamics simulation in step 2.1 is carried out in the NPT ensemble. The simulation temperature of the NPT ensemble is 300 - 310 K, the pressure is 1 atmosphere, and the particle mesh Ewald method is used to handle long-range electrostatic interactions.

[0045] The time of the conventional molecular dynamics simulation is any value within the range of 0.001 - 10 μs. Preferably, the time of the conventional molecular dynamics simulation is 0.1 - 5 μs, such as 0.001 μs, 0.1 μs, 0.5 μs, 1 μs, 1.5 μs, 2 μs, 2.5 μs, 3 μs, 3.5 μs, 4 μs, 5 μs, 6 μs, 7 μs, 8 μs, 9 μs, 10 μs.

[0046] The number of times of the conventional molecular dynamics simulation is 1 - 5 times, such as 1, 2, 3, 4, 5. Preferably, the number of times of the conventional molecular dynamics simulation is 3 times.

[0047] In a specific embodiment of the present invention, the time of the conventional molecular dynamics simulation is 1 μs, and the simulation is performed 3 times.

[0048] The heavy atom distance between the contact ligand and the binding site residue described in step 2.2 is not greater than For example, not greater than

[0049] The RMSD of the ligand and the contact probability of the binding site residue are calculated based on the heavy atoms of the amino acid residues. When the RMSD of the ligand is in a stable state after fluctuations, it can be considered that the ligand may have been in a binding state. Then, the ligand contact probability of all residues interacting with the ligand is calculated to determine the important residues. After that, the interaction map between the residues and the ligand is calculated. The starting time is the moment when the ligand interacts with the protein, and the final time is the moment when the ligand is stably bound in the pocket of the protein.

[0050] When the contact probability between the residue and the ligand is not less than 0.7, this site is considered to be the binding site of the ligand and the ligand-gated ion channel.

[0051] The method for calculating the binding free energy of the ligand and the ligand-gated ion channel described in step 2.2 includes molecular dynamics / generalized Born surface area method, thermodynamic integration method or free energy perturbation method.

[0052] The interaction between the residue and the ligand is that all heavy atoms of the residue are within a distance not greater than all heavy atoms of the ligand. When calculating the interaction map between the residue and the ligand, the interaction intensity should be normalized according to the inverse of the distance.

[0053] The starting point of the targeted molecular dynamics simulation described in step three is the representative pose structure of the ligand-gated ion channel binding to the ligand, and the end point is the crystal structure of the Holo state.

[0054] The targeted molecular dynamics simulation is carried out in the NPT ensemble.

[0055] The simulated temperature of the NPT ensemble is 300 - 310 K, and the pressure is under periodic boundary conditions of one atmosphere. Among them, for the bond lengths connected to H atoms, the limiting threshold is The cutoff value of the non-bonded interaction is

[0056] Step three mentioned above also includes performing a clustering analysis on the trajectory of the targeted molecular dynamics simulation to obtain the initial intermediate state structure.

[0057] The conventional molecular dynamics simulation described in step four is performed in NAMD2 and NAMD3.

[0058] Before performing the conventional molecular dynamics simulation in step four, it is necessary to perform energy minimization and heating on the simulation system.

[0059] The energy minimization of the simulation system includes achieving energy minimization by adding restraining forces to the ligand-gated ion channel and lipids.

[0060] The restraining force is any value within the range of 0.5 - 15 kcal.mol -1 . Preferably, the restraining force is any value within the range of 1 - 10 kcal.mol -1 . Preferably, the restraining force is any value within the range of 1 - 10 kcal.mol, such as 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 kcal.mol -1 .

[0061] The heating includes raising the temperature of the simulation system in the NVT ensemble to 300 - 350 K within a time of 100 - 300 ps. Preferably, the heating includes raising the temperature of the simulation system to 300 - 320 K within a time of 200 - 300 ps.

[0062] The step size of the conventional molecular dynamics simulation is 1 - 2 fs, and the trajectory is saved every 20 - 100 ps. Preferably, the step size of the conventional molecular dynamics simulation is 2 fs, and the trajectory is saved every 100 ps.

[0063] The conventional molecular dynamics simulation described in step four includes setting the cutoff distance of the non-bonded interaction to be any value within the range of, such as Preferably, the cutoff distance of the non-bonded interaction is

[0064] The conventional molecular dynamics simulation described in Step 4 is carried out in the NPT ensemble. The simulation temperature of the NPT ensemble is 300 - 310 K, and the pressure is 1 atmosphere. The particle mesh Ewald method is used to handle the long-range electrostatic interaction.

[0065] The time of the conventional molecular dynamics simulation described in Step 4 is any value within the range of 200 - 250 ns, such as 200 ns, 205 ns, 210 ns, 215 ns, 220 ns, 225 ns, 230 ns, 235 ns, 240 ns, 245 ns, 250 ns. Preferably, the time of the conventional molecular dynamics simulation is 210 ns.

[0066] The intermediate state stable structure described in Step 4 is the last frame of the stable trajectory of the intermediate state.

[0067] The time of the Gaussian accelerated molecular dynamics simulation described in Step 5 is any value within the range of 0.5 - 2 μs, such as 0.5 μs, 0.6 μs, 0.7 μs, 0.8 μs, 0.9 μs, 1 μs, 1.2 μs, 1.4 μs, 1.6 μs, 1.8 μs, 2 μs. Preferably, the time of the Gaussian accelerated molecular dynamics simulation is 1 μs.

[0068] Step 6 specifically includes:

[0069] Step 6.1: Integrate the intermediate state stable trajectory and the trajectory of the Gaussian accelerated molecular dynamics simulation, and extract the conformations of the ligand-gated ion channel at intervals to form a complete trajectory;

[0070] Step 6.2: Input the complete trajectory obtained in Step 6.1 into the NetworkX software to generate a dynamic network;

[0071] Step 6.3: Calculate the correlation coefficient analysis (such as the generalized correlation coefficient) between the Cα atoms of the ligand-gated ion channel residues and the ligand in the dynamic network; among them, a dynamic network is generated according to the generalized correlation coefficient and distance, and the Floyd-Warshall algorithm is used to generate an allosteric path to explain the allosteric regulation mechanism of ligand regulation of the ligand-gated ion channel;

[0072] Step 6.4: Conduct community analysis and shortest path analysis on the dynamic network.

[0073] Preferably, in the dynamic network generated in Step 6.2, the Cα atoms of the ligand-gated ion channel residues and the key atoms of the ligand are represented by nodes. If the distance between the heavy atoms of two nodes is within within, then an edge is added between the two nodes.

[0074] The thickness of each edge in the dynamic network is scaled according to the distance. The thicker the edge, the greater the correlation between the two nodes.

[0075] The path length between two nodes in the dynamic network is the sum of the individual path lengths involved between the node sets.

[0076] The residues of the ligand-gated ion channel involved in the shortest path between two nodes in the dynamic network are considered to be the residues that play an important role in the allosteric regulation of the ligand-gated ion channel.

[0077] Community analysis: The dynamic network is further divided into multiple different sub-networks using a multilevel algorithm with the generalized correlation coefficient as the edge weight. These sub-networks may represent the distribution of the communication networks of each domain of the channel, and Pymol is used for visual analysis of the community network.

[0078] Allosteric path analysis: The Floyd-Warshall algorithm is used to search for the path with the shortest distance between the allosteric site and the active site in the network. This path is often the most likely or biologically relevant signal transduction path, that is, the allosteric communication path between the functional domains of the channel protein.

[0079] The clustering analysis uses the DBSCAN algorithm and performs clustering analysis based on the backbone atoms of amino acid residues.

[0080] The ligand-gated ion channel described above includes a trimeric ligand-gated ion channel, a tetrameric ligand-gated ion channel, or a pentameric ligand-gated ion channel.

[0081] The ligand-gated ion channel described above includes an acetylcholine receptor (nAChR), a γ-aminobutyric acid receptor (GABA receptor), a glycine receptor (GlyRs), a glutamate receptor (iGluR), or a degenerin / epithelial sodium channel family (DEG / ENaC). Preferably, the degenerin / epithelial sodium channel family includes an acid-sensing ion channel (ASIC), a bile acid-sensitive ion channel (BASIC), a hydra sodium channel (HyNaC), or an FMRFamide-activated sodium channel (FaNaC).

[0082] In a specific embodiment of the present invention, the ligand-gated ion channel is an FMRFamide-activated sodium channel (FaNaC), and the ligand of the FaNaC is FMRFamide;

[0083] The binding sites of FaNaC and FMRFamide include I147, L149, I152, Y156, L168, Q172, M182, N183, I185, F188, T271, Y272, G273, V274, F453, and Y454.

[0084] The residues with a binding free energy lower than -0.5 during the binding of FMRF amide to FaNaC include I147, I152, L168, Q172, M179, M182, N183, F188, T271, and V274.

[0085] The sites where the ligand-residue contact probability is not less than 0.7 during the binding of FMRF amide to FaNaC include I147, L149, I152, and N183.

[0086] The sites where FMRF amide and FaNaC protein interact strongly after 200 ns include V274, Y272, T271, N183, and Q172.

[0087] The shortest allosteric path for the binding of FMRF amide to FaNaC is N183 - F181 - S367 - L353 - N345 - C435 - Y329, V274 - N183 - I147 - E263 - Q525 - E518, or F453 - Q455 - R461 - H209 - L203 - K197 - I185 - N191 - W387.

[0088] The key sites during the binding of FMRF amide to FaNaC include Q172, N183, V272, F453, and L563.

[0089] In the second aspect of the present invention, a molecular dynamics simulation method for the binding path and allosteric path of a ligand-gated ion channel and a ligand is provided. The molecular dynamics simulation method includes the path recognition method described in the first aspect.

[0090] In the third aspect of the present invention, a path recognition device for a ligand-gated ion channel based on molecular dynamics simulation is provided. The device includes a storage module, and the storage module includes a program and / or model for implementing the path recognition method described in the first aspect.

[0091] In the fourth aspect of the present invention, an application of the path recognition method described in the first aspect, the molecular dynamics simulation method described in the second aspect, or the path recognition device described in the third aspect in simulating the binding path and allosteric path of a ligand-gated ion channel and a ligand is provided.

[0092] In the fifth aspect of the present invention, an application of the path recognition method described in the first aspect, the molecular dynamics simulation method described in the second aspect, or the path recognition device described in the third aspect in screening allosteric regulators that bind to a ligand-gated ion channel is provided.

[0093] The allosteric regulator includes an allosteric activator or an allosteric inhibitor.

[0094] Beneficial effects:

[0095] (1) The main object of the present invention is ligand-gated ion channels. As one of the most complex membrane proteins, ion channels have characteristics such as a larger and more complex system and smaller conformational changes. Therefore, their ligand-binding path and allosteric path are also relatively complex, and they have greater reference significance for other proteins. Compared with the proteins disclosed in the prior art, the research scheme of the present invention has better scalability and extensibility.

[0096] (2) The present invention uses molecular dynamics to simulate the process of random binding of ligands to ligand-gated ion channels. Molecular dynamics simulation is the simulation method closest to experimental conditions in molecular simulation, and can simulate the microscopic evolution process of the system at the atomic level. Compared with conventional rigid docking, the method of the present invention has higher accuracy, pays attention to the conformational changes of ligand-gated ion channels during the binding process, and the obtained allosteric path is more real and effective.

[0097] (3) With the aid of molecular dynamics simulation, contact probability and ligand residue interaction analysis, free energy analysis and dynamic network analysis, the present invention identifies the ligand-binding path of LGICs, extracts the important intermediate state structures between the Apo state and the Holo state of LGICs and analyzes their allosteric paths. Taking the FMRFamide-gated sodium channel (FaNaC) as an example, the binding path of FMRFamide entering the binding pocket and its binding sites are successfully identified. The binding sites include I147, L149, I152, Y156, L168, Q172, M182, N183, I185, F188, T271, Y272, G273, V274, F453 and Y454. Compared with the prior art, more binding sites are found. By extracting the intermediate state structures between the Apo state and the Holo state, the present invention is also applicable to the research of other ligand-gated channels.

[0098] (4) The present invention uses a variety of simulation methods such as conventional molecular dynamics, targeted molecular dynamics and Gaussian accelerated molecular dynamics, which ensures the random binding process of ligands while making the extracted important intermediate state conformations more stable and credible, and overcomes the time-scale limitation of traditional molecular dynamics simulation, reducing the computing power consumption.

[0099] (5) The present invention integrates various analysis methods including dynamic network analysis to reveal allosteric pathways, intermediate state structures, and ligand binding pathways. Dynamic network analysis can analyze the transmission efficiency of allosteric information, identify allosteric pathways and their important residues that play important roles in allosteric information transmission, and assist in the research of the gating mechanism of LGICs and the development of important allosteric modulators. Taking the FMRFamide-gated sodium channel (FaNaC) as an example, starting from the identified binding sites as allosteric starting points, various possible allosteric pathways are obtained. For example, the allosteric pathway starting from N183 includes N183-F181-S367-L353-N345-C435-Y329, the allosteric pathway starting from V274 includes V274-N183-I147-E263-Q525-E518, and the allosteric pathway starting from F453 includes F453-Q455-R461-H209-L203-K197-I185-N191-W387. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] Figure 1 is a flowchart of the present invention.

[0101] Figure 2 is a flowchart of molecular dynamics simulation.

[0102] Figure 3 is the Ramachandran plot of the constructed complete protein, where A is the Apo state and B is the Holo state.

[0103] Figure 4 is the top view and side view of the simulation system.

[0104] Figure 5 is the RMSD plot of FMRFamide and FaNaC residues.

[0105] Figure 6 is the binding free energy composition plot of FMRFamide and FaNaC residues, where ΔE(Internal) is the thermodynamic energy, ΔE(Electrostatic) is the electrostatic interaction energy, ΔE(VDW) is the van der Waals energy, and ΔG(Binding) is the binding free energy of FMRFamide and FaNaC residues.

[0106] Figure 7 is the mutual contact probability plot of FMRFamide and FaNaC residues.

[0107] Figure 8 is the interaction plot of FMRFamide and FaNaC protein over time.

[0108] Figure 9 is the binding pathway plot of FMRFamide during the binding process.

[0109] Figure 10 It is the shortest allosteric pathway map of FaNaC. Specific implementation mode

[0110] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only partial embodiments of the present invention, rather than all. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the scope of protection of the present invention.

[0111] Embodiment 1

[0112] Taking the FMRFamide-activated sodium channel (FaNaC) as an example in the embodiments of this application, the conformational changes and allosteric pathway analysis methods of ligand-gated ion channels binding to ligands based on molecular dynamics simulation are described, but it is not limited to this example. This method is also applicable to the conformational changes and allosteric pathways of other ligand-gated ion channels binding to ligands. The overall analysis flowchart of this application is as Figure 1 shown, and the molecular dynamics simulation process is as Figure 2 shown.

[0113] The described method includes the following steps:

[0114] Step 1: Obtain the crystal structures of the Apo state and Holo state of FaNaC, complete the missing structures, and obtain the simulation system.

[0115] Specifically:

[0116] Obtain the Apo state structure (PDBID: 7YVC) and Holo state structure (PDBID: 7YVB) of FaNaC from the Protein Data Bank (PDB). Other components in the crystal structure except the protein are deleted from the Apo state structure, and other components in the crystal structure except the protein and the ligand of the corresponding subunit are deleted from the Holo state structure.

[0117] Use AlphaFold2 and RoseTTAFold integrated in ColabFold to predict the missing SpecificInsertionⅡ region (Apo state residue numbers: 489-508, Holo state residue numbers: 490-506), compare the pLDDT scores of the generated structures, select the structure with higher confidence, and align it to the Apo state and Holo state structures. Then use the completed structure as a template for homology modeling in SwissModel to generate an optimized structure. After the structure is generated, use the Ramachandran plot to verify the rationality of the structure, where Figure 3Through the analysis of 118 structures with a resolution above 2 Å and an r-factor not greater than 20%, in the most favorable region, the expected value of a good-quality model is above 90%. Figure 3 And Table 1 shows that the quality evaluation of the three-dimensional structure is qualified.

[0118] Table 1: Analysis Table of Ramachandran Plot Results for Apo State and Holo State

[0119]

[0120]

[0121] Construct a simulation system:

[0122] Before constructing the simulation system, add the ligand molecule FMRFamide to the Apo-state structure. A total of 3×3 ligands are placed at a distance of at least above the extracellular domain (ECD) of FaNaC and arranged randomly in a square. Then, construct the simulation system using CHARMM-GUI. The system includes the receptor in the Apo state, the newly added ligand, the phospholipid bilayer, solvent molecules, and ions. The protein structure is aligned in the Orientations of Proteins in Membranes (OPM) database and then inserted into the lipid bilayer composed of 100% phosphatidylcholine (POPC). Then, solvate it with TIP3P water molecules in the box. Add 0.15 mol / L NaCl to the aqueous phase to neutralize the system. The total number of atoms in the simulation system is approximately 300,000. The top view and side view of the simulation system are as shown in Figure 4 Figure.

[0123] Step 2: Perform conventional molecular dynamics simulations on the system after adding ligands to the Apo state of the FaNaC protein and the Holo state system.

[0124] Simulation trajectory: Perform 3×1 μs of conventional molecular dynamics simulations respectively. Then, analyze the ligand RMSD value, the contact probability between the ligand and nearby residues, and the residue-ligand interaction during the process of the ligand entering the binding pocket for each simulation to obtain the initial trajectory of the binding of the FaNaC protein to the ligand and extract the representative pose structure of the binding of the FaNaC protein to the ligand. Specifically:

[0125] cMD simulation settings:

[0126] First, add 10.0 kcal.mol -1 . and 5.0 kcal.mol -1 . The restraint force is used to minimize the energy. Then, the system is heated to 310 K in the NVT ensemble within 250 ps. After that, 10 ns of unrestrained equilibration is carried out in the NPT ensemble at 310 K and 1 bar using the Langevin thermostat and the Monte Carlo barostat. The cutoff distance for non-bonded interactions is set to The particle mesh Ewald (PME) method is used to calculate long-range electrostatic interactions. The time step is 2 fs. The CHARMM 36m force field is used for proteins, lipids, and ions, and the CHARMM TIP3P model is used for water molecules. The ligand parameters are generated using the Ligand Designer module of CHARMM-GUI. Trajectory snapshots are saved every 100 ps. All the above steps are performed in NAMD2 and NAMD3, and the initial trajectory of the binding of the FaNaC protein to the ligand is obtained through cMD simulation. The above settings are used in all cMD, tMD, and GaMD.

[0127] Initial trajectory analysis:

[0128] In the generated initial trajectory, 1 frame is extracted every 10 frames. The trajectory file is stripped of hydrogen atoms and aligned with the first frame conformation. Based on the heavy atoms of the ligand, the root mean square deviation (RMSD) is calculated using the MDTraj software package. The RMSD formula is as follows:

[0129]

[0130] N is the number of atoms calculated, and δ i is the position of the i-th atom in a certain frame minus its position in the reference conformation (position offset). The RMSD results of FMRF amide and FaNaC protein residues are as Figure 5 shown.

[0131] When the RMSD value stabilizes at a certain value after a period of time, it can be considered that the state of FMRF amide has been basically stable and may have bound to the binding pocket. Further calculate the binding free energy of the ligand to nearby residues, and use the conformation with the lowest binding free energy as the best binding pose for further analysis.

[0132] Molecular mechanics / generalized Born surface area method (MM / GBSA) is an effective tool for obtaining the binding free energy of protein-ligand interactions and protein-protein interactions. The MM / GBSA method is used to calculate the binding free energy of ligand-receptor interactions. The formula is as follows:

[0133] ΔG binding =G complex -(G receptor +G ligand ) (2)

[0134] where Gcomplex , G receptor and G ligand represent the receptor-ligand complex, the free energy of the receptor, and the free energy of the ligand, respectively. The free energy can be calculated by the following formula:

[0135] G = E gas + G sol - TS (3)

[0136] E gas = E int + E ele + E vdw (4)

[0137] G sol = G psolv + G npsolv (5)

[0138] Among them, E gas is the gas-phase energy, which is the sum of the thermodynamic energy E int , the van der Waals energy E vdw , and the electrostatic interaction energy E ele . G psolv and G npsolv are the polar and non-polar contributions to the solvation energy G sol . T is the temperature, and S is the total conformational entropy.

[0139] All free energy calculations were performed in MolAICal based on the NAMD module, and the results are as Figure 6 shown.

[0140] If the heavy atom distance between the FMRF amide and the residue is within , it can be considered that there is a contact and interaction between them. Then calculate the contact probability between the residue with at least 1 frame within the distance from the FMRF amide within . The contact probability formula is as follows:

[0141]

[0142] N is the total number of frames after the ligand RMSD stabilizes, and i is the number of frames in which the residue is in contact with the ligand. If the contact probability between the residue and the FMRF amide is not less than 0.7, it can be considered that this site may be the binding site of the FMRF amide. The contact probability between the FMRF amide and the residues of the FaNaC protein is as Figure 7 shown. Then, calculate the interaction map ( Figure 8 ) of the FMRF amide with the FaNaC protein over time starting from the moment when they interact, and extract the binding path ( Figure 9) The domains of the FaNaC protein to which the FMPF amide binds include Finger, β6-β7-loop, and Knuckle. The specific binding residues are I147, L149, I152, Y156, L168, Q172, M182, N183, I185, F188, T271, Y272, G273, V274, F453, and Y454 (see Table 2).

[0143] Table 2: Table of FMRF Amide Binding Sites

[0144]

[0145]

[0146] Cluster analysis:

[0147] Perform cluster analysis on the initial trajectory of the binding of the FaNaC protein to the ligand obtained from the cMD simulation. Specifically: Based on the RMSD of the backbone atoms of FaNaC, use the DBSCAN clustering algorithm to generate a set of MD conformations representing the cMD trajectory. All snapshots are superimposed using the protein Cα atoms to eliminate overall rotation and translation. Use the clustering function DBSCAN in sklearn for cluster analysis. Use evaluation metrics such as CHI (Calinski-Harabasz Index), SCI (Silhouette Coefficient Index), and DBI (Davies-Bouldin Index) to evaluate the clustering results. The centroid structure of each cluster in the clustering results is used as the structure of the ligand-binding representative pose of the cluster in subsequent analysis and extraction.

[0148] Step 3: Perform tMD simulation on the representative pose structure of the binding of the FaNaC protein to the FMRF amide and the Holo-state crystal structure to obtain the initial intermediate-state structure during the process from the Apo state to the Holo state.

[0149] Perform cluster analysis on the cMD trajectory in Step 2, extract the representative pose structure of the binding of the FaNaC protein to the FMRF amide and the Holo-state structure as the starting point and the ending point for tMD simulation, and perform cluster analysis on the tMD simulation trajectory again to extract the initial intermediate-state structure. Specifically:

[0150] tMD targeted molecular dynamics simulation:

[0151] After the above cMD simulation optimization, in the dynamics simulation module of the molecular dynamics simulation software, use the NPT ensemble to perform targeted molecular dynamics simulation in the optimized simulation system to obtain the tMD simulation trajectory.

[0152] Specifically, first, the G559-E577 and R91-T107 of the FaNaC transmembrane domain in two states are aligned respectively, and an external force is applied to this part of the amino acids, and tMD simulation is performed on this part of the structure. In the present invention, the tMD simulation uses the NPT ensemble and is implemented in the TMD module in NAMD. The starting and ending coordinates of the tMD simulation are the representative pose structure of the above-obtained FaNaC protein after binding a ligand in the Apo state and the crystal structure of the Holo state respectively. The simulation temperature of the NPT ensemble is 310K, the pressure is the periodic boundary condition of 1 atmosphere, the temperature control uses the Berendsen temperature coupling method, and the pressure control uses the isotropic correction method based on molecules. The SHAKE algorithm is used to limit the bond lengths of all bonds connected to H atoms, and the limiting threshold is The cutoff value of the non-bonded interaction is The PME method is used to handle the long-range electrostatic interaction. In the Extra Parameters, the added TMD elastic constant k is The principle formula of the TMD simulation is as follows:

[0153]

[0154] Among them, RMS(t) is the instantaneous best-fit RMS distance between the current coordinates and the target coordinates, and RMS * (t) linearly evolves from the initial RMSD of the first TMD step to the final RMSD of the last TMD step. The elastic constant k is scaled according to the number N of target atoms.

[0155] After that, clustering analysis is performed, and the DBSCAN clustering centroid structure is used as the representative structure to extract the initial intermediate state structure.

[0156] Step 4: Perform cMD simulation on the initial intermediate state structure generated by the tMD simulation again to obtain the stable trajectory of the intermediate state, and use the last frame as the stable structure of the intermediate state.

[0157] Perform an unconstrained equilibrium of 210 ns cMD simulation to obtain the stable trajectory of the intermediate state. Calculate the acceleration parameters according to the cMD simulation trajectory of the last 10 ns, and select the final structure of the cMD simulation as the starting structure for the subsequent GaMD simulation.

[0158] Step 5: Perform Gaussian accelerated molecular dynamics simulation on the stable structure of the intermediate state to obtain the Gaussian accelerated molecular dynamics simulation trajectory.

[0159] Specifically:

[0160] Gaussian accelerated molecular dynamics simulation:

[0161] Gaussian Accelerated Molecular Dynamics (GaMD) can enhance the conformational sampling of biomolecules by adding a harmonic gain potential to the system when the system potential V(r) is lower than the reference energy E. The specific formula is as follows:

[0162]

[0163] Among them, k is the harmonic force constant, and E and k are adjustable parameters automatically determined by three enhanced sampling principles. The reference energy E can be calculated by formula (10):

[0164]

[0165] Among them, V min and V max are the minimum and maximum potential energies of the system respectively. To ensure the effectiveness of formula (10), k must satisfy k ≤ 1 / (V max - V min ).

[0166] A total of 1 μs of GaMD simulation was performed to obtain the Gaussian accelerated molecular dynamics simulation trajectory.

[0167] Step 6: Perform dynamic network analysis on the intermediate state stable structure obtained in Step 4 and the Gaussian accelerated molecular dynamics simulation trajectory obtained in Step 5 to obtain the allosteric pathway of ligand binding to the FaNaC protein.

[0168] Specifically:

[0169] Using the dynetan and NetworkX software packages, all MD simulation trajectories (a total of 3 μs) of each system are used to generate a dynamic network for each system. In the dynamic network, the Cα atoms of the receptor residues and the key atoms of the ligand are represented by nodes. If the distance between the heavy atoms of two nodes is within within 75% of the sampling time, an edge is added between the two nodes. The generalized correlation coefficient r MI is derived from the mutual information I calculated from the positions of a pair of nodes, which defines the information transfer between two nodes i and j at a given simulation time and can be calculated by the following formula:

[0170]

[0171] Among them, N is the total number of simulation frames, is the digamma function, n i and n jis the number of frames in which the positions of nodes i and j are close to their positions in the reference (frame 0), and they are averaged by varying the reference frame across all simulations. After that, the mutual information is converted to the generalized correlation coefficient r by applying Equation (12). MI 。

[0172]

[0173] where d = 3 for the (x, y, z) dimensions that describe each node.

[0174] The calculation of the mutual information estimate between a pair of nodes can be described more specifically as follows: Given a pair of nodes i and j and a reference simulation frame f0, the positions of each node are compared with their positions in all other simulation frames. For each frame, the largest change in the x, y, and z dimensions is selected to represent the "distance" of the node from its position in frame f0 (the largest change among all dimensions is used, rather than the Cartesian distance). For each frame in the simulation, then the largest distance of either i or j is selected and used to sort all the frames. Taking the k nearest neighbors in the "simulation frame" space, i.e., the k frames in which the positions of nodes i and j have the least change compared to frame f0, the largest changes between the x, y, and z dimensions of each node i and j can be determined, giving d i and d j 。These two distances are used as cut-off values to select the frames in which nodes i and i are closer to their respective positions in frame f0 than d i and d j 。n i and n j are the number of simulation frames that meet these criteria. From the first frame to the last frame, the same calculation is performed with f0 varying, resulting in the average calculated in Equation (11), and then the generalized correlation coefficient is achieved by applying the mutual information estimate in Equation (12).

[0175] The thickness of each edge in the dynamic network is scaled by distance, and thicker edges indicate greater correlation. Then the original network is further divided into different sub-networks, called communities, using the louvain algorithm. Nodes within a community have stronger connections compared to nodes in other communities.

[0176] The shortest path uses the Floyd-Warshall algorithm to search for the path with the shortest distance between two nodes. The shortest path is often the most probable or biologically relevant path, thus revealing potential allosteric regulatory paths. The path length between two nodes in the dynamic network is equal to the sum of the individual path lengths involved between the sets of nodes. The shortest allosteric path of the FaNaC protein is as Figure 10 shown.

[0177] The effect of ligands on the communication of various regions of the FaNaC protein can be revealed by the connection strength of each community in community analysis. The length of the shortest path can reflect the communication strength between two distal regions of the FaNaC protein, and the path length is inversely proportional to the communication strength. Moreover, the residues involved in the shortest path can be considered as residues that play an important role in the allosteric regulation of the FaNaC protein. The allosteric path lengths in different states are shown in Table 3, where bound1 refers to the ligand randomly binding to 1 subunit, bound2 refers to the ligand randomly binding to 2 subunits, and bound3 refers to the ligand randomly binding to 3 subunits.

[0178] Table 3: Allosteric path lengths in different states

[0179] Apo Holo bound1 bound2 bound3 Q172(I) → F120(I) 51.27865 56.90098 48.92580 52.49364 54.08831 N183(I) → Y329(I) 52.94625 63.09850 55.34131 58.97180 57.68794 V274(I) → E512(I) 51.79076 54.78561 52.41488 61.08676 61.21520 F453(I) → W387(I) 58.84540 60.68897 56.18272 57.53765 58.72909 F453(III) → Y329(I) 34.72402 31.43227 26.10906 27.71235 38.15702 Q172(II) → F120(II) 62.36751 59.15418 43.82961 48.05642 53.53212 N183(II) → Y329(II) 67.67898 64.92993 47.73903 53.72912 56.93255 V274(II) → E512(II) 59.18195 60.31535 48.20278 57.70681 68.78790 F453(II) → W387(II) 59.52714 70.47295 50.91729 57.00602 62.56134 F453(I) → Y329(II) 37.54015 30.55141 44.77317 41.08413 30.31419 Q172(III) → F120(III) 63.13228 67.89224 52.24861 54.16419 47.36536 N183(III) → Y329(III) 63.20590 61.24998 50.15932 56.35562 54.81283 V274(III) → E512(III) 58.71634 69.59874 51.59504 53.68820 55.49185 F453(III) → W387(III) 68.80867 55.75729 49.11062 50.12463 54.92790 F453(II) → Y329(III) 48.19439 46.23229 28.95703 34.78877 33.40814

[0180] In summary, this embodiment achieves the following:

[0181] (1) Identifying the ligand-binding path and ligand-binding sites of FaNaC;

[0182] (2) Extracting the important intermediate state structures between the Apo state and the Holo state of FaNaC;

[0183] (3) Revealing the allosteric path between ligand binding and the pore domain of FaNaC.

[0184] Comparative Example 1

[0185] Compared with Example 1, step 2 was removed and step 3 was directly performed using tMD simulation. The results were not applicable to ligand-gated ion channels because there are multiple transition states between the Apo state and the Holo state, and the time scale they span is very long. Moreover, the system of ligand-gated ion channel proteins is large and complex, and direct use of tMD simulation cannot accurately capture the ligand-binding process and conformational changes.

[0186] Therefore, in Example 1, first using cMD simulation to allow the ligand to autonomously bind to the ion channel is not only for studying the ligand-binding process, but also for making the channel tend to change towards the Holo state, reducing the difficulty of subsequent simulation in crossing the energy barrier.

[0187] Comparative Example 2

[0188] Replacing the analysis method in step 2 of Example 1 with a single RMSD, contact probability, or binding free energy to determine the ligand-binding state, the results show that the credibility of determining whether the ligand binds to the channel protein through a single index is low.

[0189] Therefore, in Example 1, determining whether the ligand binds to the channel protein through three indicators of RMSD, contact probability, and binding free energy is more credible than a single indicator.

[0190] Comparative Example 3

[0191] Compared with Example 1, Step 3 was removed, and after 2 cMD simulations, GaMD simulation was directly performed. The result was that the conformational changes of the ligand-gated ion channel could not be captured. This was because there were still multiple transition states between the obtained resting-state ligand-bound structure and the Holo state structure. Due to the large size and complex structure of the ion channel system, it was difficult to cross such a long time scale only using enhanced sampling algorithms.

[0192] Therefore, in Example 1, in order to study the allosteric pathway of LGICs, tMD simulation was used to apply external forces between some amino acids in a single domain of the protein to target the protein from one state to another, so as to capture the conformational changes from the ligand-binding site of the finger domain to the distal transmembrane domain.

[0193] Comparative Example 4

[0194] Compared with Example 1, Step 4 was removed, and after tMD simulation, GaMD simulation was directly performed. The result was that the stability of the ligand-gated ion channel decreased during a longer simulation time.

[0195] Therefore, performing cMD simulation again after tMD simulation in Example 1 could provide a longer time scale to observe the dynamic behavior of the system, including the vibration, conformational changes of the ion channel protein, and solvent effects, etc.; moreover, cMD could also verify whether the structure obtained by tMD simulation remained stable during a longer simulation time, thereby verifying the stability of the structure.

[0196] Among them, 2 cMDs were used in Example 1, which were performed before and after tMD simulation respectively, effectively improving the stability of the protein structure, bringing it to a relaxed state, and enhancing the stability and credibility of the structure. The first cMD was used after pre-equilibration, aiming to release the channel from the constrained state and allow the ligand to bind to the ligand-gated channel in a as natural a manner as possible; the second cMD was used after tMD simulation, aiming to further stabilize the structure after tMD simulation and make the structure in a relaxed state.

[0197] Comparative Example 5

[0198] Compared with Example 1, Step 5 was removed, and GaMD was not performed after the second cMD simulation. The result showed that although the intermediate state structures of some domains had been sampled by tMD simulation, within a limited simulation time, the system could not fully explore all its possible conformational spaces.

[0199] Therefore, in Example 1, after the second cMD simulation, GaMD simulation is carried out. By adding a harmonic boosting potential to smooth the biomolecular potential energy surface and reduce the energy barrier, compared with other enhanced sampling algorithms, its advantage lies in that there is no need to set predefined reaction coordinates or collective variables. GaMD simulation provides unconstrained enhanced sampling, jumping out of the energy barrier to explore the conformational space as much as possible. It is beneficial to simulate complex biological processes and is very suitable for complex membrane proteins such as ligand-gated ion channels. The above is only the specific implementation manner of the present invention, but the protection scope of the present invention is not limited to this embodiment. Other modifications and implementation manners designed by those skilled in the art based on this invention should also be covered within the protection scope of the present invention.

Claims

1. A path identification method for ligand-gated ion channels based on molecular dynamics simulation, characterized in that: The method includes: Step 1: Obtain the crystal structures of the Holo state and Apo state of the ligand-gated ion channel as the initial structure and build a simulation system; Step 2: Perform conventional molecular dynamics simulation on the simulation system to obtain the initial trajectory of the ligand-gated ion channel binding to the ligand and extract the representative posture structure of the ligand-gated ion channel binding to the ligand; Step 3: Perform targeted molecular dynamics simulation on the representative posture structure of the ligand-gated ion channel and the crystal structure of the Holo state to extract the initial intermediate state structure; Step 4: Perform conventional molecular dynamics simulation on the initial intermediate state structure to obtain the intermediate state stable trajectory and extract the intermediate state stable structure; Step 5: Perform Gaussian accelerated molecular dynamics simulation on the intermediate stable structure to obtain the Gaussian accelerated molecular dynamics simulation trajectory; Step 6: Perform dynamic network analysis on the intermediate state stable trajectory obtained in step 4 and the Gaussian accelerated molecular dynamics simulation trajectory obtained in step 5 to obtain the allosteric pathway of the binding of the ligand to the ligand-gated ion channel.

2. The path identification method according to claim 1, characterized in that: The simulation system includes the crystal structure of the Apo state or the crystal structure of the Holo state of the ligand-gated ion channel; the simulation system also includes a ligand, a phospholipid bilayer, solvent molecules and ions.

3. The path identification method according to claim 1, characterized in that: The distance between the ligand-gated ion channel and the ligand is above.

4. The path identification method according to claim 1, characterized in that: The distance between the extracellular domain of a ligand-gated ion channel and the ligand is above.

5. The path identification method according to claim 2, characterized in that: The total number of atoms in the simulation system is 200,000-400,000.

6. The path identification method according to claim 1, characterized in that: The step 2 specifically includes: Step 2.1: Minimize the energy and heat the simulation system, and perform unconstrained conventional molecular dynamics simulation under the NPT ensemble; Step 2.2: Calculate the root mean square deviation of the ligand heavy atoms, the contact probability between the ligand and the binding site residues, and the binding free energy between the ligand and the ligand-gated ion channel to determine the ligand binding site; Step 2.3: Perform cluster analysis on the initial trajectory of the ligand-gated ion channel binding to the ligand, and extract the centroid structure of each cluster as the representative posture structure of the ligand-gated ion channel binding to the ligand.

7. The path identification method according to claim 6, characterized in that: The heating described in step 2.1 involves raising the temperature of the simulated system in the NVT ensemble to 300-350 K within a period of 100-300 ps.

8. The path identification method according to claim 7, characterized in that: The heating includes raising the temperature of the simulated system to 300-320K within 200-300ps.

9. The path identification method according to claim 6, characterized in that: The time of conventional molecular dynamics simulation in step 2.1 is 0.001-10 μs.

10. The path identification method according to claim 9, characterized in that: The duration of the conventional molecular dynamics simulation in step 2.1 is 1 μs.

11. The path identification method according to claim 6, characterized in that: The number of conventional molecular dynamics simulations in step 2.1 is 1-5 times.

12. The path identification method according to claim 11, characterized in that: The number of conventional molecular dynamics simulations in step 2.1 is 3.

13. The path identification method according to claim 1, characterized in that: The targeted molecular dynamics simulation described in step three starts with a representative posture structure of the ligand-gated ion channel bound to the ligand, and ends with a crystal structure of the Holo state.

14. The path identification method according to claim 1, characterized in that: The intermediate state stable structure described in step 4 is the last frame of the stable trajectory of the intermediate state.

15. The path identification method according to claim 1, characterized in that: The time of conventional molecular dynamics simulation in step 4 is 200-250ns.

16. The path identification method according to claim 15, characterized in that: The time of conventional molecular dynamics simulation in step 4 is 210 ns.

17. The path identification method according to claim 1, characterized in that: The number of conventional molecular dynamics simulations in step 4 is 1.

18. The path identification method according to claim 1, characterized in that: The step six comprises: Step 6.1: Integrate the intermediate state stable trajectory and the trajectory of Gaussian accelerated molecular dynamics simulation, extract the conformation of the ligand-gated ion channel at intervals, and form a complete trajectory; Step 6.2: Input the complete trajectory obtained in step 6.1 into the NetworkX software to generate a dynamic network; Step 6.3: Calculate the correlation coefficient between the Cα atoms of the ligand-gated ion channel residues in the dynamic network and the ligand; Step 6.4: Perform community analysis and shortest path analysis on the dynamic network.

19. The path identification method according to claim 18, characterized in that: The correlation coefficient analysis is a generalized correlation coefficient.

20. The path identification method according to claim 1, characterized in that: The ligand-gated ion channels include acetylcholine receptors (nAChR), gamma-aminobutyric acid receptors (GABA receptors), glycine receptors (GlyRs), glutamate receptors (iGluRs) or degenerative protein / epithelial sodium channel family (DEG / ENaC).

21. The path identification method according to claim 20, characterized in that: The degenerate protein / epithelial sodium channel family includes acid-sensing ion channel (ASIC), bile acid-sensing ion channel (BASIC), hydra sodium channel (HyNaC) or FMRF amide-activated sodium channel (FaNaC).

22. The path identification method according to claim 1, characterized in that: The ligand-gated ion channel is FaNaC, the ligand is FMRF amide, and the sites at which FMRF amide binds to FaNaC include one or more of I147, L149, I152, Y156, L168, Q172, M182, N183, I185, F188, T271, Y272, G273, V274, F453 or Y454.

23. The path identification method according to claim 22, characterized in that: The shortest allosteric pathways for FMRF amide binding to FaNaC include N183-F181-S367-L353-N345-C435-Y329, V274-N183-I147-E263-Q525-E518, or F453-Q455-R461-H209-L203-K197-I185-N191-W387.

24. A path identification device based on molecular dynamics simulation of ligand-gated ion channels, characterized in that: The device includes a storage module, and the storage module includes a program and / or model for implementing the path identification method described in any one of claims 1-23.

Images