A molecular dynamics simulation method for enhanced sampling

Abstract

This invention discloses an enhanced sampling molecular dynamics simulation method, belonging to the field of molecular simulation technology. The specific steps include: S1, defining reaction coordinates; S2, transitioning from the binding state to the encounter state; S3, determining constraints on the conformation of the protein and / or ligand; S4, transitioning from the encounter state to the dissociation state; and S5, result analysis. When using umbrella sampling to obtain the protein-ligand binding free energy, this invention proposes dividing the protein-ligand binding process into two parts: from the dissociation state to the encounter state, and then from the encounter state to the binding state. It employs universal reaction coordinates, does not impose constraints on the relative positions of the protein and ligand, and uses multiple sequentially performed umbrella sampling to sample the high-probability binding paths of the protein and ligand. This invention provides a molecular dynamics method that can quickly, accurately, and with high automation calculate the binding free energy of flexible proteins and flexible ligands, as well as the actual binding / dissociation paths.

Description

Technical Field

[0001] This invention relates to the field of molecular simulation technology, and in particular to a molecular dynamics simulation method with enhanced sampling. Background Technology

[0002] In the field of biomolecular interaction research, theoretical computational methods have become as important as experimental methods. Among theoretical computational methods, molecular dynamics simulations are a crucial tool for studying the mechanisms of action of biomolecular molecules at the atomic level. However, given the complexity of biomolecular molecules, this method has limitations in the size of the systems it can study and the length of the observation time. To obtain a more comprehensive record (i.e., sampling) of the behavior of biomolecular molecules within a limited observation time (i.e., simulation time), it is usually necessary to use enhanced sampling methods to accelerate traditional molecular dynamics simulations. Commonly used methods include umbrella sampling and metadynamics, which have been widely used in the study of protein conformational changes and the binding processes of proteins and small molecule ligands.

[0003] The basic idea of ​​umbrella sampling is to divide the system's conformation space into a series of windows along one or more response coordinates and sample in each window. In each window, a bias potential energy (also known as umbrella energy) is applied to restrict the sampling of the system conformation to that window, thereby increasing the sampling efficiency of system conformations that have a low probability of occurring in specific regions on the response coordinates.

[0004] However, the existing methods have the following drawbacks: (1) Due to the restriction on the orientation of the ligand relative to the protein, the ligand can only bind / dissociate with the protein in a specific orientation, and the obtained binding / dissociation path often deviates from the real path; (2) In order to correct the binding free energy error caused by the restriction, six additional umbrella sampling simulations are required, which makes the process too complicated and not conducive to the promotion of the method; (3) The reaction coordinates used need to be manually defined according to the specific research system, and there are many human-computer interaction operations; (4) If there are other important dimensions that are not considered in the simulation, and you want to add restrictions to that dimension, you need to re-perform the umbrella sampling simulation, which leads to high trial and error costs; (5) Since the selected reaction coordinate RMSD has ambiguity in characterizing the conformation space (conformations with very different characteristics often have the same RMSD value), it is easy to cause a large error in the average potential (PMF) along the RMSD obtained by umbrella sampling; (6) The conformation space of the protein-ligand binding state lacks a clear and strict definition, which may lead to errors in obtaining the related calculations of the binding state.

[0005] To address the aforementioned shortcomings, this invention aims to develop a molecular dynamics method capable of rapidly, accurately, and with high automation calculating the binding free energy of flexible proteins and flexible ligands, as well as the actual binding / dissociation pathways. Summary of the Invention

[0006] The present invention aims to provide a molecular dynamics simulation method with enhanced sampling to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] A molecular dynamics simulation method with enhanced sampling, comprising the following steps:

[0009] S1. Define reaction coordinates: Perform a conventional molecular dynamics simulation of medium length on the protein-ligand complex, and determine the amino acids at the binding interface based on the obtained trajectory, i.e., the amino acids whose distance to the ligand is less than a specific value. Select the distance between the centroid of the amino acid at the binding interface and the ligand and / or the native contact of the protein-ligand as the reaction coordinates.

[0010] S2, from bound state to encountered state: Using the natural contact as the reaction coordinate, run multiple umbrella-shaped samplings of medium length. Stop sampling when the ligand is detected to have largely left the protein.

[0011] S3. Determine the constraints between protein and / or ligand conformations: Analyze the path of the ligand leaving the protein in the previous step to determine the structure of the main complex encountered. Here, we take the ligand as an example for explanation. Protein conformation constraints are similar.

[0012] S4. Encounter state to dissociation state: Based on the RMSD constraints added to the protein and ligand respectively obtained in the previous step, multiple medium-length umbrella-shaped samplings are run using centroid distance as the reaction coordinate;

[0013] S5. Results Analysis: Based on the umbrella sampling simulation of the protein ligand binding process, combined with the simulation results of protein-complex and protein or ligand individual systems, a multidimensional free energy surface is constructed through principal component analysis and statistical physics principles to obtain the protein ligand binding free energy and binding / dissociation pathway.

[0014] Preferably, in step S3, the ligand structure can be restricted to the vicinity of its conformation where it binds to the protein. The threshold of the restriction parameter is determined by analyzing the conformational changes of the ligand in the previous step of the medium-length simulation. The ligand dissociation trajectory is projected onto the two-dimensional free energy surface of the ligand obtained in the second part. The position of the ligand conformation in the encountered complex is determined by comparing different dissociation paths, thereby determining the range of conformational restriction on the ligand fragment.

[0015] Preferably, the reaction coordinates mentioned in step S1 also include DRMSD, i.e., RMSD based on the distance between amino acids and / or bases.

[0016] Preferably, sampling of the protein-ligand binding process can be replaced by enhanced sampling schemes that can basically guarantee the integrity of the binding pathway, such as stretching dynamics, meta-dynamics, and weighted ensemble.

[0017] Preferably, in order to reduce the required sampling conformation space, the funnel potential energy used in Funnel Metadynamics can be used to further reduce the required sampling space, while ensuring that the pathway of ligand-binding proteins is not severely disturbed.

[0018] Compared with the prior art, the present invention has the following beneficial effects:

[0019] This invention proposes dividing the protein-ligand binding process into two parts: from the dissociated state to the encountered state, and then from the encountered state to the bound state. From the encountered state to the bound state (or equivalently from the bound state to the encountered state), the resulting binding / dissociation pathway more closely approximates the actual pathway because no additional constraints are added to the system. Simultaneously, due to the interaction between the protein and ligand, the conformational changes of the protein and ligand from the bound state to the encountered state are relatively small, requiring less sampling space, and ensuring sufficient sampling of the conformational space even without restrictions. While the conformations of the protein and ligand are restricted separately during the process from the encountered state to the dissociated state, no restrictions are added to their relative positions. Therefore, the relative orientation of the protein and ligand can change, maximizing the fidelity of the binding / dissociation pathway.

[0020] Since there are multiple possible paths for binding / dissociation, this invention uses multiple sequentially performed umbrella-shaped sampling to sample the protein-ligand binding path. Only the path with the highest probability (not necessarily the highest probability or all possible paths) can be selected for sampling, thereby reducing the required sampling range.

[0021] This invention uses common reaction coordinates such as RMSD, centroid distance, and native contact. These reaction coordinates can be automatically generated using Bash scripts and are applicable to other protein-ligand systems, thereby greatly reducing manual operation and improving the level of automation.

[0022] Building upon existing simulations, new constraints can be easily added to enhance sampling of potentially undersampled dimensions, and then a more accurate free energy surface can be constructed using statistical physics principles. This can accelerate computation and avoid wasting computational resources.

[0023] Principal component analysis provides the principal components as coordinates to characterize the conformation space, which increases the discriminative power of the conformations and provides a more accurate description of the conformation space. Attached Figure Description

[0024] Figure 1A schematic diagram of a novel thermodynamic cycle for calculating the protein-ligand binding free energy (ΔG°bind) with added constraints;

[0025] Figure 2 This is a schematic diagram of the sampling method;

[0026] Figure 3 This is a schematic diagram showing the distribution of the path obtained from umbrella sampling from the binding state to the encounter state on the first two principal components of the PCA.

[0027] Figure 4 This is a schematic diagram of the free energy surfaces on the first two principal components of the PCA obtained from umbrella-shaped sampling from the binding state to the encounter state.

[0028] Figure 5 This is a schematic diagram of the free energy surfaces on the first two principal components of the PCA obtained from umbrella-shaped sampling from the encounter state to the dissociation state;

[0029] Figure 6 This is a schematic diagram of the free energy surfaces of a protein on its two principal components in the PCA head.

[0030] Figure 7 This is a flowchart of a molecular dynamics simulation method that enhances sampling. Detailed Implementation

[0031] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments:

[0032] The specific implementation process is as follows:

[0033] like Figure 1As shown, this technique also employs a thermodynamic cycling method with added constraints to calculate the protein-ligand binding free energy, but it differs significantly from existing techniques. First, the binding process of the protein and ligand is divided into two parts: the protein and ligand first come into contact to form an encounter complex or an encounter state, and then the binding of the protein and ligand becomes increasingly tight, eventually forming a stable binding state. Second, when using umbrella sampling to calculate the binding free energy of the protein and ligand under constrained conditions, no additional constraints are imposed on the relative positions of the protein and ligand (except for the constraints added to each window of the umbrella sampling). The calculation of each free energy change is as follows: Without additional constraints, umbrella sampling is used to calculate the free energy (ΔGb) of the protein and ligand from the encounter state to the bound state; constraints (corresponding to rp and rL) are added to the RMSD of the protein and ligand respectively, and then umbrella sampling is used to calculate the free energy (ΔGrE) of the protein and ligand from the dissociated state to the encounter state under the constraints; the complete conformational space is collected by performing traditional long-term molecular dynamics simulations on the protein and ligand systems separately (ligands such as RNA oligomers generally have small molecular weights and no fixed structure). Long-term simulations can obtain a more accurate conformational space; sampling of the protein conformational space can also be based on umbrella sampling of the RMSD of all protein backbone atoms, and then the free energy changes ΔGrP and ΔGrL caused by the added constraints are calculated based on statistical physics principles (i.e., two molecular dynamics simulations); among them, ΔGr' can be calculated based on the relevant sampling results of the umbrella sampling simulation for calculating ΔGrE and ΔGb and combined with statistical physics principles, while ΔGv can be directly calculated by analytical formulas. Therefore, the calculation of these two free energy changes does not require additional simulations.

[0034] A novel thermodynamic cycle for calculating the protein-ligand binding free energy (ΔG°bind) by adding constraints. rp and rL represent the constraints added to the protein and ligand, respectively. ΔGrP and ΔGrL represent the changes in free energy of the dissociated system after adding constraints rp and rL, respectively; ΔGrE represents the change in free energy of the protein and ligand from the dissociated state to the encounter state under the presence of constraints; ΔGr' represents the change in free energy of the protein-ligand encounter state after adding constraints rp and rL; ΔGb represents the change in free energy of the protein and ligand from the encounter state to the bound state; and ΔGv represents the change in free energy caused by the change in system concentration due to placing the protein and ligand in a specific space.

[0035] The following describes the specific implementation scheme (i.e., calculation steps) for calculating the free energy (ΔGrE and ΔGb) of the two parts of the protein ligand binding free energy using umbrella sampling, as follows: Figure 7 As shown:

[0036] (1) Define reaction coordinates: Perform a moderate-length conventional molecular dynamics simulation of the protein-ligand complex. Based on the obtained trajectory, determine the amino acids at the binding interface, i.e., amino acids whose distance to the ligand is less than a specific value (e.g., 0.6 Å). Select the distance between the centroid of the amino acid at the binding interface and the ligand, and / or the native contacts of the protein and ligand, as the reaction coordinates. Automatically generate the definition file of the reaction coordinates for the simulation software (here, GROMACS software and PLUMED plugin) using a Bash script (NC.sh).

[0037] (2) From bound state to encountered state: Multiple umbrella samplings of moderate length were run using the natural contact as the reaction coordinate. Sampling was stopped when the ligand was detected to have largely left the protein. Multiple umbrella samplings were used to enhance the collection efficiency of the binding / dissociation pathway. Here, umbrella sampling was performed sequentially window by window, starting from the window closest to the bound state and moving towards the encountered state one by one; after completing the previous window, the process moved to the next window.

[0038] Figure 3 The distribution of the umbrella-shaped sampling paths from the binding state to the encounter state on the two principal components of the PCA head. The trajectories can be divided into two categories. Seven trajectories pass through the lower main channel, and three trajectories (i.e., parallel simulations 2, 3, and 7) pass through from above. The encounter state is identified as the location with the most dissociated trajectory conformations visited, located at (-2.5, -0.5).

[0039] Figure 4 The free energy surfaces on the two principal components of the PCA head, obtained by umbrella-shaped sampling from the bound state to the encounter state. The encounter state is located at (-2.5, -0.5) with a free energy of 27.0 kJ / mol; the bound state is located at (2.7, 1.8) with a free energy of -15.1 kJ / mol.

[0040] (3) Determine the constraints on protein and / or ligand conformations: Analyze the ligand exit path from the protein in the previous step to determine the structure of the main encountering complex; this is explained using ligands as an example, and protein conformation constraints are similar. Because ligands are flexible, they undergo significant conformational changes after dissociation from the protein, making complete sampling of the ligand conformational space in a protein-ligand complex system impractical. Assume that the conformational change occurring the instant the ligand leaves the protein-binding interface is relatively small, while most conformational changes occur after dissociation. Therefore, the ligand structure can be constrained to the vicinity of its protein-binding conformation. The threshold of constraint parameters (such as RMSD) is determined by analyzing the ligand conformational changes in the medium-length simulation in the previous step. Specifically, the ligand dissociation trajectory is projected onto the two-dimensional free energy surface of the ligand obtained in the second part, and the position of the ligand conformation (i.e., conformation I) in the encountering complex is determined by comparing different dissociation paths. See [link to explanation]. Figure 2This allows for the determination of the range of conformational constraints on ligand fragments. A Bash script is used to automatically generate the necessary configuration file for these constraints based on the umbrella sampling trajectory.

[0041] Figure 2 The diagram illustrates a method to accelerate sampling by narrowing the sampling space of ligand conformations through restrictions without affecting the ligand binding process. The gray lines represent the hypothetical two-dimensional free energy surface of the ligand in aqueous solution. A and B are two stable conformations of the ligand (for simplicity, A is assumed to be the most stable conformation). C is the conformation when the ligand binds to the protein to form a complex, and I is the conformation at the moment the ligand immediately separates from the protein (i.e., encounters the complex). The black arrow from C to I indicates the path of conformational change from C to I near the protein-protein binding interface; the black arrow from I to A indicates the path of conformational change from I to A in aqueous solution; and the red arrow from C to A indicates the path of conformational change from C to A in aqueous solution. Due to the presence of the protein, the path of conformational change from C to A may differ from the path without the protein. By restricting the ligand conformation to the blue dashed box, molecular dynamics simulations can still capture the dynamic path of the ligand dissociating from the protein without needing to capture ligand conformations outside the blue dashed box. Furthermore, because the ligand conformation is confined to the vicinity of conformation I, umbrella sampling simulations can relatively easily capture the pathway where the ligand rebinds to the protein to form an encounter complex. Without these constraints, most of the simulation time would likely be spent sampling the two stable conformations A and B, making it difficult to obtain conformation I, where the ligand forms the encounter complex near the protein interface.

[0042] (4) Encounter state to dissociation state: Based on the RMSD constraints added to the protein and ligand respectively obtained in the previous step (completed using the Bash script RMSD.sh), multiple medium-length umbrella samplings are run using the centroid distance (completed using the Bash script COM.sh) as the reaction coordinates.

[0043] Figure 5 The free energy surface on the first two principal components of the PCA, obtained from umbrella-shaped sampling from the encounter state to the dissociation state. The encounter state is located at (-2.3, -2.0), with a free energy of 16.4 kJ / mol; the dissociation state is widely distributed, appearing as a plane with almost equal energy on the potential energy surface, with a free energy of 50.0 kJ / mol. Note that the first two principal components here are... Figure 3 and Figure 4 The differences are as follows;

[0044] Figure 6 : The free energy surface of the protein on the two principal components of its PCA head. The projections of the average conformation of the protein in the bound state (origin in the figure) and the encountered state (× in the figure) onto this free energy surface are marked.

[0045] (5) Results analysis: Based on the umbrella sampling simulation of the protein ligand binding process, combined with the simulation results of protein-complex and protein or ligand individual systems, a multidimensional free energy surface is constructed through principal component analysis and statistical physics principles to obtain the protein ligand binding free energy and binding / dissociation pathway.

[0046] The following example, using the RNA methylation recognition protein YTHDC1 and the RNA oligomer system, demonstrates the application results of the enhanced sampling method proposed in this invention:

[0047] Table 1 lists all molecular dynamics simulation types and durations, with a total of 5 simulations and a total duration of 7.5 microseconds. Based on simulation trajectory 4, and using protein-ligand conformational characteristics, the encounter state was determined using a Bash script. The resulting dissociation paths and encounter states are shown in [Table 1]. Figure 3 The free energy surface constructed through simulation results and principal component analysis (PCA) is shown below. Figure 4-6 . Figure 1 The calculation results of the changes in various free energies required to calculate the binding free energy are shown in Table 2.

[0048] Table 1: List of all simulations recognizing the YTHDC1 protein and RNA oligomer system (ns represents nanoseconds)

[0049]

[0050] Table 2: Figure 1 Calculation results of the changes in various free energies required to calculate the combined free energy

[0051] Energy Item Energy value (kcal / mol) <![CDATA[ΔG r P ]]> 4.0 <![CDATA[ΔG r L ]]> 12.2 <![CDATA[ΔG v ]]> 20.7 <![CDATA[ΔG b ]]> -42.1 <![CDATA[ΔG r E ]]> -35.6 <![CDATA[ΔG r’ ]]> 3.7

[0052] The above descriptions are merely embodiments of the present invention, and common knowledge such as specific technical solutions and / or characteristics are not described in detail here. It should be noted that those skilled in the art can make various modifications and improvements without departing from the technical solutions of the present invention, and these should also be considered within the scope of protection of the present invention. These modifications and improvements will not affect the effectiveness of the implementation of the present invention or the practicality of the patent. The scope of protection claimed in this application should be determined by the content of its claims, and the specific embodiments described in the specification can be used to interpret the content of the claims.

Claims

1. A molecular dynamics simulation method with enhanced sampling, characterized in that, The specific steps include: S1. Define reaction coordinates: Perform a conventional molecular dynamics simulation of medium length on the protein-ligand complex, and determine the amino acids at the binding interface based on the obtained trajectory. When the distance from the amino acid to the ligand is less than a specific value, select the distance between the centroid of the amino acid and the ligand at the binding interface and / or the native contact of the protein-ligand as the reaction coordinates. S2, from bound state to encountered state: Using the natural contact as the reaction coordinate, run multiple umbrella-shaped samplings of medium length. Stop sampling when the ligand is detected to have largely left the protein. S3. Determine the constraints between protein and / or ligand conformations: Analyze the pathways by which ligands leave the protein in the previous step to determine the structures of the main complexes encountered. In step S3, the ligand structure is restricted to the vicinity of the conformation where the ligand binds to the protein. The threshold of RMSD restriction is determined by analyzing the conformational changes of the ligand in the medium-length simulation in the previous step. The ligand dissociation trajectory is projected onto the two-dimensional free energy surface of the ligand obtained in step S3. The position of the ligand conformation in the encountered complex is determined by comparing different dissociation paths, thereby determining the range of conformational restriction of the ligand fragment. S4. Encounter state to dissociation state: Based on the RMSD constraints added to the protein and ligand respectively obtained in the previous step, multiple medium-length umbrella-shaped samplings are run using centroid distance as the reaction coordinate; S5. Results Analysis: Based on the umbrella sampling simulation of the protein ligand binding process, combined with the simulation results of protein-complex and protein or ligand individual systems, a multidimensional free energy surface is constructed through principal component analysis and statistical physics principles to obtain the protein ligand binding free energy and binding / dissociation pathway.

2. The molecular dynamics simulation method with enhanced sampling according to claim 1, characterized in that: The reaction coordinates mentioned in step S1 also include RMSD based on the distance between amino acids and / or bases.

3. The molecular dynamics simulation method with enhanced sampling according to claim 1, characterized in that: The sampling of the protein-ligand binding process is replaced with enhanced sampling schemes that ensure the integrity of the binding pathway, such as stretching dynamics, meta-dynamics, and weighted ensemble.

4. The molecular dynamics simulation method with enhanced sampling according to claim 1, characterized in that: To reduce the required sampling conformation space, the funnel potential used by FunnelMetadynamics is employed to further reduce the required sampling space while ensuring that the pathway of ligand-binding proteins is not severely disturbed.

Images