A method for identifying key residues in small molecule-protein binding

Abstract

The present invention discloses a method for identifying key residues for exploring small molecule-protein binding based on stretching molecular dynamics. The present invention simulates a system composed of different small molecules and proteins by stretching molecular dynamics method, studies the dissociation kinetics of small molecules, and analyzes the key factors affecting the dissociation rate of small molecules. The method can study mutants of viruses and analyze the inhibitory effects of inhibitors on different viral mutations by means of computational biology simulation. Compared with traditional molecular dynamics methods, this method can accelerate the simulation to the physiological reaction time scale. The method can provide sufficient theoretical basis for structure-based drug design and improve the inhibitory effect of inhibitors on viral mutants from the atomic level.

Description

Technical Field

[0001] The present invention relates to the field of computational drug design using molecular dynamics methods, and in particular to a method for identifying key residues of small molecule-protein binding based on stretched molecular dynamics. Background Art

[0002] Inhibitors have significant inhibitory effects on proteases that have clear targets in some viruses and are related to their physiological processes. However, viruses have obvious mutability, which can cause proteases to produce drug-resistant mutations to inhibitors, which makes it difficult to develop effective inhibitors that can deal with drug-resistant mutants. The current research and development cycle of inhibitors is long, and the long-term use of inhibitors will have great limitations. In order to improve the effectiveness of inhibitors and enhance their ability to face mutation risks, the focus of drug design is on the binding and dissociation kinetics of the interaction between inhibitors and proteases. Studies have shown that most inhibitors currently exhibit differential kinetics for different proteases, and the potential inhibition mechanism in this process is crucial for protein structure-based drug design. Different inhibitors have different binding rates and dissociation rates for different proteases, and exploring the factors that affect this rate constant is an important driving force for drug design.

[0003] The effectiveness of an inhibitor is usually measured using the equilibrium constant K d This value can be expressed by the dissociation rate k off and the binding rate k on The lower the equilibrium constant, the better the inhibitory effect. Therefore, in the design process of inhibitors, how to increase the binding rate and reduce the dissociation rate is the key goal in the design process. Molecular dynamics (MD) simulation is an effective tool for studying the dynamic process of inhibitors and proteins. Molecular dynamics is based on Newton's second law and can solve the motion process of a large number of atoms and molecules in the simulation system. This method is widely used in biology, computational chemistry, and pharmacy. In drug design, this method can identify the key factors that affect the binding and dissociation process of inhibitors in the structure-based drug design process at the atomic and molecular levels. Studies have shown that the size, topological structure, and stiffness of inhibitors can affect the binding process of different inhibitors. Compared with the binding process, people have studied the dissociation process relatively less. Mutations in proteases can often reduce the effectiveness of inhibitors by significantly increasing the dissociation rate. Therefore, mutations in proteases have a significant effect on the dissociation rate, and this mechanism is of great significance for avoiding inhibitor failure.

[0004] All-atom simulation can effectively track the atomic details of the dynamics of the simulated system. Studying the dissociation process of inhibitors through all-atom simulation methods is of great significance for optimizing the interactions in the inhibitor design process. The time scale of the dissociation process of inhibitors under physiological conditions reaches the order of milliseconds and above, and the time scale of classical all-atom molecular dynamics simulation is usually in the order of nanoseconds to microseconds. Therefore, the simulation of the dissociation process of inhibitors under physiological conditions exceeds the limit of classical molecular dynamics simulation. Steered molecular dynamics (SMD) simulation can apply a force to the system, so steered molecular dynamics simulation can accelerate the dissociation process of inhibitors. This method has been widely used to study systems such as ligand-receptor interactions, protein-protein interactions, and protein unfolding processes. In addition, this method can also calculate the free energy diagram of intermolecular interactions. Through the stretched molecular dynamics method, the potential mechanism that affects the dissociation rate of inhibitors can be found more efficiently in atomic details, which can effectively improve the research and development process of new inhibitors. In summary, introducing computational biology methods into the field of drug design can assist drug developers in designing targeted drugs and improve the tolerance of inhibitors to viral mutation behavior. Summary of the invention

[0005] Based on the above shortcomings of the prior art, the technical problem solved by the present invention is to propose a method for identifying key residues in small molecule-protein binding based on stretching molecular dynamics. The present invention uses the stretching molecular dynamics method as a simulation method to simulate the small molecule-protein system, analyzes the number of hydrogen bonds and interatomic distances of the simulated trajectory, analyzes the key interaction residues in the wild type; analyzes the effect of the wild type on the inhibitor effect after mutation; uses umbrella sampling to calculate the binding energy and quantitatively analyzes the binding strength of small molecules and proteins. This method can assist in the design of targeted drugs and improve the tolerance of inhibitors to viral mutations.

[0006] In order to solve the above technical problems, the present invention is implemented by the following technical solutions:

[0007] The identification of key residues for small molecule-protein binding based on stretched molecular dynamics includes the following steps:

[0008] S1: Obtain protein crystal structures from the PDB database and fill in missing crystal structures; download small molecule structure files from the small molecule database; generate small molecule force field files and topology files; and construct a complex simulation system;

[0009] S2: Energy minimization, annealing simulation, pre-equilibration under NVT and NPT ensembles are performed on the simulation system, and the constraints are released to perform a shorter conventional molecular dynamics simulation to generate the initial structure of stretched molecular dynamics;

[0010] S3: Perform pre-equilibrium under the NPT ensemble, remove constraints at different stretching rates to perform stretching molecular dynamics; select reaction coordinates for umbrella sampling; extract key indicators of small molecule-protein interaction and calculate binding free energy.

[0011] Wherein, S1 comprises the following steps:

[0012] S11: Download the appropriate protein crystal structure from the PDB database and select the crystal structure with complete structure and high resolution;

[0013] S12: Download the molecular structure file from the small molecule database and optimize the structure;

[0014] S13: Generate small molecule force field files and topology files using CGenFF tool;

[0015] S14: Construction of small molecule-protein complex system.

[0016] Furthermore, S2 includes the following steps:

[0017] S21: Use the CHARMM36 force field to characterize the interactions between protein molecules and atoms, and use the TIP3P water molecule model to solvate the entire system;

[0018] S22: Add ions to the solvated system to make the net charge of the system zero; Na + and Cl - Neutralize the system charge; and set the osmotic pressure to 150 mM;

[0019] S23: perform annealing simulation on the simulated system, heating the system from 3K to 310K within 100ps;

[0020] S24: Use the V-rescale algorithm to couple the temperature of the simulation system and perform a 100ps simulation under the NVT ensemble;

[0021] S25: The pressure of the simulation system is coupled using the Parrinello-Rahman algorithm, and a 100 ps simulation is performed under the NPT ensemble;

[0022] S26: A conventional molecular dynamics simulation is performed on the pre-equilibrated system after being released from restraint for about 50 ns.

[0023] Finally, S3 includes the following steps:

[0024] S31: Extraction of simulated structures for stretched molecular dynamics from conventional molecular dynamics;

[0025] S32: build a simulation box according to the system size;

[0026] S33: The Parrinello-Rahman algorithm is used to couple the pressure of the simulation system and a 100 ps simulation is performed under the NPT ensemble;

[0027] S34: Generate configuration files with springs of different stiffness coefficients and extension rates;

[0028] S35: Derivation to obtain the expression of pulling force;

[0029] S36: Perform stretching molecular dynamics simulation until the small molecule falls off from the binding pocket;

[0030] S37: Extracting the initial conformation of the umbrella sampling simulation from the stretched molecular dynamics simulation trajectory;

[0031] S38: sampling at a certain distance along the reaction coordinate;

[0032] S39: Using simple harmonic potential to limit the distance to the normal line, simulate each window until convergence;

[0033] S310: Calculate the free energy curve by weighted histogram.

[0034] Therefore, the present invention has the following advantages:

[0035] 1. This paper sets the research target as the main protease of coronavirus, constructs a small molecule-protein complex system, places the small molecule near the binding site during the simulation, applies a simple harmonic potential to the center of mass of the small molecule, observes the kinetic process of small molecule dissociation, provides an effective target for drug design, and can greatly improve the drug design cycle.

[0036] 2. Compared with traditional molecular dynamics methods, tensile molecular dynamics can accelerate the dissociation process by applying a force to small molecules. In addition, computational biology can be used to achieve multi-site and multi-type mutations in proteins, and simulation methods can be used to study targets that can limit and enhance the inhibitory effects of small molecules. It can greatly improve the research on tolerance to viral mutations and provide multi-faceted assistance for drug development. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 is a flow chart of the method of the present invention;

[0038] Figure 2 It is a schematic diagram of a simulation system of the present invention;

[0039] Figure 3 The tension curves and distance curves along the reaction coordinate of different small molecules in the present invention;

[0040] Figure 4 The curves of hydrogen bond numbers and interatomic distances of different small molecules in the present invention;

[0041] Figure 5 is the free energy of small molecules in different systems of the present invention. DETAILED DESCRIPTION

[0042] In order to explain the technical solution of the present invention in more detail, the technical solution of the present invention is specifically described below in conjunction with embodiments and drawings.

[0043] See also Figure 1 , a method for identifying key residues in small molecule-protein binding based on stretched molecular dynamics, including the following steps:

[0044] S1: Download the crystal structure of the protein from the PDB database, obtain the small molecule structure file, and construct the small molecule-protein complex structure;

[0045] S2: Establish a conventional molecular dynamics model, and the simulation software uses Gromacs-2018.3 version;

[0046] S3: Establish stretch exchange molecular dynamics simulation. The simulation software uses Gromacs-2018.3 version, the control software plumed version is 2.8.0, and the spring stiffness coefficient is set to 100 kJ / mol·nm 2 The stretching rate was set to 0.1nm / ps to 10nm / ps, with a total of 3 orders of magnitude. Umbrella sampling was performed according to the reaction coordinates and the average force potential curve was plotted.

[0047] Further, S1 includes the following steps:

[0048] S11: Download the crystal structure of the protein from the PDB database and delete the original ligand in the crystal structure;

[0049] S12: Download the small molecule structure from the small molecule database and dock the small molecule to the location of the original ligand;

[0050] S13: Using molecular modeling software to perform homology modeling on the monomers to generate a symmetrical dimer structure;

[0051] S14: Generate small molecule force field parameters and topology files using CGenFF tool;

[0052] S15: Construction of small molecule-protein complex structures.

[0053] Further, S2 includes the following steps:

[0054] S21: Use the pdb2gmx command to generate the structural topology file and position constraint file. The force field is selected according to the software version corresponding to the CHARMM36 force field. In this simulation, CHARMM36 Mar-2019 is selected, and the solvent is the TIP3P water model.

[0055] S22: Generate a simulation box using the editconf command and place the protein structure in the center of the box;

[0056] S23: Use the solvate command to add water solvent to the box;

[0057] S24: Use the genion command to adjust the system to electrical neutrality and set the ion concentration;

[0058] S25: Minimize the energy of the system by adjusting

[0059] emtol=10.0

[0060] nsteps=500000

[0061] Set the step size of the energy minimization process and the maximum force on the atoms of the system when convergence;

[0062] S26: Perform annealing simulation on the system, by adjusting

[0063] annealing_time=0 100 0 100

[0064] annealing_temp=3 310 3 310

[0065] Set the system's initial temperature and heating temperature, and set the time required for this heating process;

[0066] S27: pre-balance the system, including two rounds of pre-balancing under NVT and NPT ensembles;

[0067] S271: Pre-balance the system under the NVT ensemble, in the configuration file through

[0068] tcoupl=V-rescale

[0069] Set up the system to couple the temperature using the V-rescale algorithm;

[0070] S272: Pre-balance the system under the NVT ensemble, in the configuration file through

[0071] pcoupl=berendsen

[0072] Set up the system to couple the pressure using the Berendsen algorithm;

[0073] S28: After removing the constraints, perform a regular molecular dynamics simulation, in the configuration file

[0074] integrator=md

[0075] nsteps=25000000

[0076] It shows that after the system is freed from constraints, the leapfrog algorithm is used for integration, and the simulation time is set to 50ns.

[0077] Furthermore, S3 includes the following steps:

[0078] S31: extract the structure with a certain relaxation time from the conventional molecular dynamics simulation by trjconv command as the initial structure of the stretched molecular dynamics;

[0079] S32: pre-equilibrate the complex for 100 ps in the NPT ensemble;

[0080] S33: Generate simulation profiles with different stiffness coefficients and stretch rates through scripting

[0081]

[0082]

[0083] S34: During the stretching process, a constant pulling force acts on the center of mass of the small molecule

[0084] In the configuration file, two groups are set in a reaction coordinate system, and simple harmonic waves are used to stretch the virtual atoms at the center of mass of the small molecule along the set direction;

[0085] pull_ncoords=1

[0086] pull_ngroups=2

[0087] pull_group1_name=UNL

[0088] pull_group2_name=Protein

[0089] pull_coord1_type=umbrella

[0090] pull_coord1_geometry=direction

[0091] S35: The expression of this force is f=k spring (vt-v0), where k springrepresents the spring stiffness coefficient, v represents the stretching rate, and t is the simulation time; taking the derivative of this formula, we can get f′=k spring v;

[0092] The rate of change of the tension can be calculated to span multiple orders of magnitude;

[0093] f′ min =0.002nm / ps*200kJ / mol·nm 2 =0.67pN / nm

[0094] f′ max =0.02nm / ps*2000kJ / mol·nm 2 =67pN / nm

[0095] S36: Perform stretching molecular dynamics simulation until the small molecule falls off from the active site;

[0096] gmx mdrun-deffnm pull.tpr-pullx pullx-pullf pullf

[0097] S37: Repeat the simulation under the same conditions for multiple times;

[0098] S38: drawing a tensile rate-fracture force curve according to the simulation trajectory;

[0099] S39: Find the moment when the curve has a significant drop from the rupture force curve, and analyze the residues that affect the small molecule-protein binding strength based on the trajectory;

[0100] S310: extracting an initial configuration for umbrella sampling simulation from a stretching molecular dynamics simulation trajectory;

[0101] S311: sampling at certain intervals, with the reaction coordinates set as the distance between the small molecule and the center of mass of the binding site;

[0102] echo 0|gmx_mpi trjconv-s pull.tpr-f pull.xtc-o conf.gro-sep

[0103]

[0104] S3111: Calculate the reaction coordinates frame by frame through the script. In this case, the sample results are as follows: 0 0.505 1 0.506 2 0.509 3 0.503 4 0.501 5 0.505 6 0.505 7 0.510 8 0.501 9 0.509 10 0.507

[0115] … 490 5.266 491 5.264 492 5.293 493 5.288 494 5.289 495 5.338 496 5.366 497 5.355 498 5.349 499 5.329 500 5.352

[0127] S3112: Set the window size of the sampling distance, usually set to 0.1-0.2nm;

[0128] S3113: extracting the configuration that meets the window size setting, and performing the extraction;

[0129] S312: performing a short-time simulation on the extracted configuration under the NPT ensemble;

[0130] S313: limiting the distance to the normal line through the simple harmonic potential, and simulating each window for a certain time until the PMF curve converges;

[0131] S314: Plotting free energy curves by weighted histograms

[0132] The fracture force curve simulated by this method can better reflect the residues that have a greater contribution to the dissociation kinetics of small molecules. By anchoring the interaction of small molecules on the main chain of residues, the tolerance of inhibitors to viral mutations can be improved to a certain extent. Therefore, this implementation of the enhanced inhibitory effect on viral mutation based on stretching molecular dynamics can better assist medical personnel in structure-based drug design.

[0133] The above are only specific implementation methods of the present invention, but the protection scope of the present invention is not limited thereto, and any changes or substitutions that are not conceived through creative work should be covered within the protection scope of the present invention.

Claims

1. A method for identifying key residues in small molecule-protein binding, comprising the following steps: S1: Download the protein crystal structure co-crystallized with the inhibitor small molecule from the PDB database to fill in the missing crystal structure; download the inhibitor small molecule structure file with existing structure from the small molecule database; generate the small molecule force field file and topology file; and construct the complex simulation system; S2: Minimize the energy of the simulation system, heat the system from 3K to 310K within 100ps for annealing simulation, perform 100ps pre-equilibration in NVT and NPT ensembles respectively, release the constraints and perform 50ns of conventional molecular dynamics simulation to generate the initial structure of stretched molecular dynamics; S3: The pressure of the simulation system is coupled using the Parrinello-Rahman algorithm, and pre-balance is performed under the NPT ensemble. The spring stiffness coefficient is set to 100 kJ / mol·nm 2 The stretching rate was set to 0.1nm / ps to 10nm / ps, a total of 3 orders of magnitude, to simulate the system, and to release the constraints to perform stretching molecular dynamics; the reaction coordinates were selected for umbrella sampling; the key indicators of the interaction between small molecules and proteins were extracted and the binding free energy was calculated.

2. The method for identifying key residues for exploring small molecule-protein binding according to claim 1, characterized in that: S1 includes the following steps: S11: Download appropriate protein crystal structures from the PDB database and select protein crystal structures that have co-crystallization with inhibitor small molecules; S12: Download the inhibitor molecular structure file with existing structure from PDB database and optimize the structure; S13: Generate small molecule force field files and topology files using CGenFF tool; S14: Construction of small molecule-protein complex system.

3. The method for identifying key residues in small molecule-protein binding according to claim 1, wherein S2 The following steps are involved: S21: Use the CHARMM36 force field to characterize the interactions between protein molecules and atoms, and use the TIP3P water molecule model to solvate the entire system; S22: Add ions to the solvated system to make the net charge of the system zero; Na + and Cl - Neutralize the system charge; and set the osmotic pressure to 150 mM; S23: perform annealing simulation on the simulated system, heating the system from 3K to 310K within 100ps; S24: Use the V-rescale algorithm to couple the temperature of the simulation system and perform a 100ps simulation under the NVT ensemble; S25: The pressure of the simulation system is coupled using the Parrinello-Rahman algorithm, and a 100 ps simulation is performed under the NPT ensemble; S26: A conventional molecular dynamics simulation is performed on the pre-equilibrated system after being released from restraint for about 50 ns.

4. The method for identifying key residues in small molecule-protein binding according to claim 1, wherein S3 The following steps are involved: S31: Extraction of simulated structures for stretched molecular dynamics from conventional molecular dynamics; S32: build a simulation box according to the system size; S33: The Parrinello-Rahman algorithm is used to couple the pressure of the simulation system and a 100 ps simulation is performed under the NPT ensemble; S34: Spring stiffness coefficient is set to 100 kJ / mol·nm 2 , the stretching rate is set to 0.1nm / ps to 10nm / ps, a total of 3 orders of magnitude, and the system is simulated. After the simulation is realized, the stretching rate is reduced and the stiffness coefficient is adjusted. The simulation is repeated to generate a configuration file; S35: Derivation to obtain the expression of pulling force; S36: Perform stretching molecular dynamics simulation until the small molecule falls off from the binding pocket; S37: Extracting the initial conformation of the umbrella sampling simulation from the stretched molecular dynamics simulation trajectory; S38: sampling at a certain distance along the reaction coordinate; S39: Using simple harmonic potential to limit the distance to the normal line, simulate each window until convergence; S310: Calculate the free energy curve by weighted histogram.

5. The method for identifying key residues for exploring small molecule-protein binding according to claim 4, characterized in that: In S31, the sudden change in the tension curve during the stretching process indicates that the interaction between the small molecule and the protein has changed, and the key interacting residues can be obtained by combining the trajectory.

6. The method for identifying key residues for exploring small molecule-protein binding according to claim 4, characterized in that: In S37, sufficient initial structures need to be extracted from the stretching molecular dynamics simulation for umbrella sampling, and the umbrella sampling effect is determined by the intersection of the weighted histogram. In S38, the reaction coordinates should be able to reflect significant changes in the system morphology. For small molecule-protein complex systems, the distance between the small molecule and the binding site can be used as the reaction coordinate, and the window size is set to 0.2nm.

Images