Method for generating base-specific nucleic acid molecule force field parameters using a reweighting algorithm

Abstract

The present application relates to a method for generating base-specific nucleic acid molecule force field parameters based on a reweighting algorithm, whether the four-nucleotide conformation generated by molecular dynamics simulation is a false intercalation conformation, and the non-bond parameters of nucleic acid nucleoside heavy atoms are specifically adjusted and introduced into the grid-based energy correction (CMAP) parameters. The present application is characterized in that, in order to solve the problem of insufficient accuracy of the existing RNA molecule force field, the reweighting algorithm is used to specifically adjust the non-bond parameters of the nucleoside heavy atoms and introduce the CMAP energy term to generate base-specific nucleic acid molecule force field parameters, improve the simulation accuracy of the force field, and test multiple RNA systems through long-time simulation, replica exchange (REMD), and solute tempering replica exchange (REST2), to further verify the effectiveness of the method and the force field parameters.

Description

TECHNICAL FIELD

[0001] The present application relates to a method for generating base-specific nucleic acid molecule force field parameters by a weighting algorithm, in particular to a method for generating base-specific nucleic acid molecule force field parameters based on a re-weighting algorithm. BACKGROUND

[0002] It is very difficult to analyze RNA structure and describe its conformational assembly by traditional experimental methods. Molecular dynamics simulation (MD) has become an important means of RNA research because it can continuously sample the conformation at the atomic scale. However, there are still two major limitations of current MD simulation: first, although the development of supercomputers has continuously improved computing power, the maximum time range of current MD simulation is still limited to microseconds, which is far from enough to reproduce many biological processes (sampling problem); second, the existing energy function (force field) is not accurate enough, which may cause bias in conformation and energy sampling (force field problem).

[0003] The current RNA molecule force field has the following two major defects: first, the stability of base stacking is overestimated; second, error intercalation conformations are easily produced when simulating four-nucleotide systems. To solve the above problems, adjusting the non-bonding parameters of the atoms related to base stacking and introducing the lattice-based energy correction term (CMAP) parameter become feasible methods for optimizing the RNA molecule force field. The present application generates base-specific nucleic acid molecule force field parameters by adjusting the non-bonding parameters of the heavy atoms of nucleosides specifically based on a re-weighting algorithm and introducing the CMAP energy term. SUMMARY

[0004] The present application provides a method for generating base-specific nucleic acid molecule force field parameters based on a re-weighting algorithm, which adjusts the non-bonding parameters of the heavy atoms of nucleosides specifically and introduces the CMAP energy term to generate base-specific nucleic acid molecule force field parameters, thereby improving the simulation accuracy of the force field and further verifying the effectiveness of the method and the force field parameters by testing multiple RNA systems through long-time simulation, replica exchange (REMD), and solute tempering replica exchange (REST2).

[0005] The technical problems solved by the present application are realized by the following technical solutions:

[0006] The present application provides a method for generating base-specific nucleic acid molecule force field parameters based on a re-weighting algorithm, which adjusts the non-bonding parameters of the heavy atoms of nucleosides specifically and introduces the CMAP energy term to generate base-specific nucleic acid molecule force field parameters, thereby improving the simulation accuracy of the force field and further verifying the effectiveness of the method and the force field parameters by testing multiple RNA systems through long-time simulation, replica exchange (REMD), and solute tempering replica exchange (REST2).

[0007] Further, the method comprises the following steps:

[0008] (1) Initial molecular dynamics simulation: REST2 simulation with solute tempering replica exchange is performed for four tetranucleotide systems (AAAA, CCCC, GGGG, UUUU) based on the most commonly used nucleic acid force field ff99bsc0χOL3 parameters;

[0009] (2) Trajectory analysis: all conformations of the simulation reference replica trajectory are extracted, the base stacking order is analyzed using DSSR software, and a script is written to distinguish whether it is a false intercalation conformation;

[0010] (3) Re-weighting adjustment of non-bonding parameters: based on whether the conformation is an intercalation structure, the non-bonding parameters of the heavy atoms of the nucleosides on the four bases are adjusted specifically;

[0011] (4) Second molecular dynamics simulation: REST2 simulation with solute tempering replica exchange is performed for five tetranucleotide systems (AAAA, CCAAU, CCCC, GACC, UUUU) based on the adjusted parameters;

[0012] (5) Trajectory analysis: all conformations of the simulation reference replica trajectory are extracted, the stacking is analyzed using DSSR software, a script is written to distinguish whether it is a false intercalation conformation, and the non-bonding energy between the heavy atoms of the nucleosides is calculated;

[0013] (6) Re-weighting adjustment of CMAP parameters: based on whether the conformation is an intercalation structure, the ζ / α backbone dihedral angle grid energy term (CMAP) is introduced by re-weighting algorithm;

[0014] (7) Obtain a set of base-specific accurate molecular force field parameters (BSFF1) and carry out testing.

[0015] Further, the self-compiled script is used to analyze the DSSR software calculation results to judge whether it is a false intercalation conformation.

[0016] Further, the re-weighting algorithm is used to adjust the non-bonding parameters of the heavy atoms of the nucleosides based on the proportion of false intercalation conformations.

[0017] Further, the re-weighting algorithm is used to introduce the ζ / α backbone grid energy correction parameter (CMAP) based on the proportion of false intercalation conformations to obtain the final nucleic acid base-specific molecular force field parameters (BSFF1).

[0018] The present application has the beneficial effects that:

[0019] (1) The application provides a method for generating base-specific nucleic acid molecule force field parameters through a reweighting algorithm, in particular, a method for generating base-specific nucleic acid molecule force field parameters based on a reweighting algorithm, aiming at the problem of insufficient accuracy of existing RNA molecule force fields, the method specifically adjusts the non-bonding parameters of nucleotide heavy atoms through a reweighting algorithm and introduces CMAP energy terms to generate base-specific nucleic acid molecule force field parameters, thereby improving the simulation accuracy of the force field, and the method is further verified to be effective through long-time simulation, replica exchange (REMD), and solute tempering replica exchange (REST2) for testing multiple RNA systems.

[0020] (2) The embodiment of the application fully considers the structural property differences of various bases and can effectively improve the simulation effect of nucleic acid molecules.

[0021] (3) The embodiment of the application adopts a method for generating nucleic acid base-specific molecule force field parameters through a reweighting algorithm, which belongs to the field of methods, and specifically designs to reweight the ff99bsc0χOL3 force field simulation of four nucleotides to determine whether the conformation is an intercalated structure, specifically adjusts the non-bonding parameters of different base nucleotide heavy atoms, and then simulates multiple four-nucleotide systems using the adjusted parameters, and finally introduces the ζ / α backbone dihedral CMAP energy term through the reweighting algorithm to achieve the purpose of generating base-specific RNA molecule force field parameters. BRIEF DESCRIPTION OF DRAWINGS

[0022] In order to more clearly illustrate the technical solutions in the embodiments of the application or the prior art, the drawings needed in the following embodiment or prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the application, and other drawings can be obtained by those skilled in the art without creative labor.

[0023] Figure 1 A method for generating base-specific nucleic acid molecule force field parameters based on a reweighting algorithm according to an embodiment of the application is shown in the flowchart.

[0024] Figure 2 A four-nucleotide system test result graph in an embodiment of the application is shown in the graph.

[0025] Figure 3 A riboswitch system test result graph in an embodiment of the application is shown in the graph.

[0026] Figure 4A schematic diagram of a stem-loop hairpin system in an embodiment of the application for ab initio folding;

[0027] Figure 5 A schematic diagram of a stem-loop hairpin system in an embodiment of the application for ab initio folding. DETAILED DESCRIPTION

[0028] In order to make the technical means, creative features, purposes and effects of the present application easy to understand, the present application will be further described below in combination with specific drawings.

[0029] The embodiment of the present application aims at the problem of insufficient accuracy of the existing RNA molecular force field, adjusts the non-bonding parameters of the nucleotide heavy atoms specifically through a reweighting algorithm, introduces a CMAP energy item to generate base-specific nucleic acid molecular force field parameters, improves the simulation accuracy of the force field, and tests multiple RNA systems through long-time simulation, replica exchange (REMD), and solute tempering replica exchange (REST2) to further verify the effectiveness of the method and the force field parameters.

[0030] The present application provides a method for generating base-specific nucleic acid molecular force field parameters based on a reweighting algorithm, which adjusts the non-bonding parameters of the nucleic acid nucleotide heavy atoms specifically and introduces grid-based energy correction (CMAP) parameters based on whether the four-nucleotide conformation generated by molecular dynamics simulation is a false intercalated conformation.

[0031] The embodiment of the present application provides a method for generating base-specific molecular force field parameters (BSFF1) of nucleic acids based on a reweighting algorithm. The method flow is divided into two steps: first, the four-nucleotide molecular dynamics simulation trajectory is preprocessed to distinguish whether the simulated conformation is a false intercalated conformation, and the non-bonding parameters of each base nucleotide heavy atom are specifically adjusted based on this using a reweighting method; second, four-nucleotide molecular dynamics simulation is performed on the adjusted parameters to distinguish whether the simulated conformation is a false intercalated conformation, and grid-based energy correction (CMAP) parameters for optimizing the main chain dihedral angle ζ / α are introduced based on this using a reweighting method, to finally generate base-specific molecular force field parameters of nucleic acids.

[0032] The embodiment of the present application adjusts the non-bonding parameters of the nucleic acid nucleotide heavy atoms and introduces CMAP parameters for the main chain dihedral angle ζ / α through a reweighting method, to reduce the generation of four-nucleotide false intercalated conformation, improve the accuracy of molecular force field simulation, and test multiple nucleic acid molecules through common simulation, replica exchange (REME), and solute tempering replica exchange method (REST2) and other enhanced sampling methods, to further verify the effectiveness of the method.

[0033] The method for generating base-specific nucleic acid molecular force field parameters through the above-mentioned reweighting algorithm specifically includes the following steps:

[0034] (1) Initial molecular dynamics simulation

[0035] To improve sampling efficiency and ensure sampling accuracy, we used the enhanced sampling method of replica exchange with solute tempering (REST2) to obtain as comprehensive a conformational space as possible with limited computing resources. REST2 simulations were performed for the AAAA, CCCC, GGGG, and UUUU four-nucleotide systems using GROMACS 2018.8 and plumed 2.6.2 software. The preprocessing of the simulation system is as follows: first, for the standard conformation of the four-nucleotide structure, the pdb2gmx tool provided by the Gromacs software was used to generate the topology file, and Tip3p water model and resistant ions were added to simulate the initial structure of the natural random protein under physiological conditions. Water molecules were added to simulate the electrically neutral solvent environment, and the steepest descent method was used for energy minimization of the system for 4000 steps. Then the system was heated to 275K for 5000 steps of NVT and 25000 steps of NPT. During the entire process, all covalent bonds connected to hydrogen atoms were constrained using the LINCS algorithm, and the Particle Mesh Ewald (PME) algorithm was used to calculate the long-range electrostatic interactions. The cutoff values for Lennard-Jones interactions and electrostatic interactions were set to (1nm). After the preprocessing, replica exchange with solute tempering (REST2) simulation was started, with 16 exchange replicas spanning a temperature range of 275K to 500K. According to the principle of REST2, the temperature of different replicas i is set according to formula (1):

[0036]

[0037] where T i represents the temperature of replica i, T0 is the lowest temperature, which is also the temperature of the reference replica, i.e. replica 0, T max represents the highest temperature, and n represents the number of replicas. According to the above formula, the temperatures of the 16 replicas are calculated as 275.00K, 286.19K, 297.82K, 309.93K, 322.53K, 335.64K, 349.29K, 363.49K, 378.27K, 393.65K, 409.66K, 426.32K, 443.65K, 461.69K, 480.46K, and 500.00K. The simulation time for a single replica is set to 200ns, with a simulation step of 2fs. A frame of results is output every 50ps, and a replica exchange is attempted every 1000 steps, i.e. 2ps.

[0038] (2) Trajectory analysis and feature calculation

[0039] For the initial molecular dynamics simulation trajectories, the trajectory of the reference copy (i.e. the copy corresponding to the lowest temperature) was extracted, the simulation structure was extracted frame by frame, and the base stacking order was analyzed using the DSSR software and a self-compiled script was used to determine whether it was an incorrect intercalation conformation. At the same time, the distance between each pair of nucleotide heavy atoms was calculated, and the non-bonding interaction energy was converted.

[0040] (3) Re-weighting adjustment of non-bonding parameters

[0041] The re-weighting method is based on the converged conformational ensemble under the existing force field parameters. After adjusting the force field parameters, each conformation is re-weighted to calculate the proportion of each conformation in the new ensemble and thus obtain its average properties without re-simulation. The calculation formula is as follows:

[0042]

[0043] where < > represents the average value, A λ′ and A λ represent the ensemble properties obtained by molecular dynamics simulation under force field parameters λ' and λ, respectively, E λ′ and E λ represent the potential energy of each conformation calculated under parameters λ' and λ, respectively, and β is a constant, which is the inverse of the product of the Boltzmann constant and the thermodynamic temperature.

[0044] Using re-weighting, the force field parameters can be linked to simulation properties that are not explicitly related to the parameters. In this way, the simulation properties under the new parameters can be estimated from the properties obtained under the current parameters, and the reasonableness of the new parameters can be evaluated. In this work, we used the re-weighting method to correct the nucleotide heavy atom non-bonding parameters for the four four-nucleotide systems based on whether each frame structure was an intercalation conformation and the corresponding non-bonding interaction energy between nucleotide heavy atoms. Four sets of base-specific nucleotide heavy atom non-bonding parameters were generated accordingly.

[0045] (4) Secondary molecular dynamics simulation

[0046] Based on the force field parameters obtained by correcting the base-specific nucleotide heavy atom non-bonding parameters in step 3, we performed secondary molecular dynamics simulations on the five four-nucleotide systems: AAAA, CCCC, UUUU, CAAU, and GACC. The simulation used the same REST2 method as step 1, and the simulation time and condition settings were the same as step 1.

[0047] (5) Trajectory analysis and feature calculation

[0048] For the molecular dynamics simulation trajectory obtained in step 4, the trajectory of the reference copy (i.e. the copy corresponding to the lowest temperature) is extracted, the simulation structure is extracted frame by frame, the base stacking order thereof is analyzed using the DSSR software, and whether it is a wrong intercalation conformation is judged by using a self-compiled script. Meanwhile, the ζ / α value of each frame structure is calculated, and the distribution thereof is counted.

[0049] (6) Adjustment of ζ / α main chain lattice energy correction parameter (CMAP) by reweighting

[0050] Based on whether each frame structure obtained by calculation and analysis in step 5 is a wrong intercalation and the distribution of the ζ / α value thereof, the method of reweighting mentioned in step 3 is used to correct the ζ / α main chain lattice energy correction parameter (CMAP), a set of general CMAP parameters is generated, and thus the final nucleic acid base-specific molecular force field parameter (BSFF1) is obtained, and the overall development process is as shown in Figure 1

[0051] The method for generating a base-specific nucleic acid molecular force field parameter based on a reweighting algorithm provided in the embodiments of the present application is used to obtain the final nucleic acid base-specific molecular force field parameter BSFF1.

[0052] Referring to Figure 1 The method for generating a base-specific nucleic acid molecular force field parameter based on a reweighting algorithm provided in the embodiments of the present application is used to obtain the final nucleic acid base-specific molecular force field parameter BSFF1.

[0053] (7) Molecular dynamics simulation test

[0054] The REST2 simulation method mentioned in step 1 is used to test the simulation effect of the BSFF1 parameter in a four-nucleotide system and compare it with the ff99bsc0χOL3 force field. The simulation time and conditions are consistent with those in step 1. The conventional molecular dynamics simulation method is used, the simulation time is 1 microsecond (μs), and the remaining simulation conditions and pretreatment methods are consistent with those in step 1. The simulation effect of the ff99bsc0χOL3 force field and the BSFF1 force field in a riboswitch system is tested. The replica exchange (REMD) enhanced sampling method is used, the BSFF1 force field is used for de novo folding simulation of a hairpin stem loop system, the number of copies is 6, the temperatures are 275K, 284.34K, 294K, 303.99K, 314.32K and 325K respectively, and the simulation time is 500 nanoseconds.

[0055] (8) Analysis after molecular dynamics simulation

[0056] ​After the simulation, we used the trjconv module of Gromacs to export the trajectory as pdb structure, the cpptraj module of AMBER to calculate the root mean square deviation (RMSD) and the DBSCAN module to cluster the structures.

[0057] Example verification

[0058] We evaluated the accuracy of BSFF1 force field by molecular dynamics simulation of tetranucleotide and riboswitch system, and further verified the global rationality of the parameters by the ability of BSFF1 force field to fold the hairpin loop system from scratch.

[0059] 4.1 Tetranucleotide

[0060] Tetranucleotide has the advantages of small simulation calculation and rich experimental values, and is often used as a standard system to test the performance of RNA molecular force field. Research has found that existing RNA force fields are prone to produce incorrect intercalation conformations when simulating tetranucleotides, which seriously affects the accuracy of the simulation. In this work, we used the TIP3P water model, and used the most commonly used nucleic acid force field ff99bsc0χOL3 and BSFF1 force field parameters to perform REST2 enhanced sampling simulation on five tetranucleotide systems: AAAA, CCCC, UUUU, CAAU, and GACC, and cluster analysis was performed on the structures. The results are shown in Figure 2 It can be found that the ff99bsc0χOL3 force field produces a large number of incorrect intercalation structures, while this problem is significantly improved in the simulation of BSFF1, which successfully reduces the frequency of incorrect intercalation conformations and increases the frequency of reasonable conformations.

[0061] Figure 2 For the test results of tetranucleotide system, by analyzing the base stacking order of each cluster structure, it can be found that the ff99bsc0χOL3 force field produces a large number of incorrect intercalation conformations, while BSFF1 produces more reasonable conformations.

[0062] 4.2 Riboswitch

[0063] Riboswitch is a class of non-coding RNA elements located on messenger RNA (mRNA) that can bind small molecule ligands to specifically regulate gene expression and play an important role in many physiological processes. The internal structure of riboswitch is complex, and its structure maintenance and stability involve multiple interactions such as base stacking and base pairing. The simulation of riboswitch can fully evaluate the properties of each parameter of the molecular force field. We performed 1 microsecond (μs) of long-time conventional molecular dynamics simulation on the 5kh8 riboswitch system under the ff99bsc0χOL3 force field and BSFF1 force field, and performed structure clustering and comparison with experimental structures, the results are shown inFigure 3 As shown, the BSFF1 force field can maintain its stable experimental structure, and the clustering results are better than those of the ff99bsc0χOL3 force field.

[0064] Figure 3 The graph shows the test results of the ribo-switching system. Figure 3 In the diagram, B represents the ff99bsc0χOL3 force field, and C represents the BSFF1 force field. By performing structural clustering and comparing it with the experimental structure, it was found that the BSFF1 force field can maintain the experimental structure, and its deviation from the experimental results (RMSD) is smaller than that of the ff99bsc0χOL3 force field.

[0065] 4.3 The stem-ring hairpin system is folded from the top.

[0066] RNA stem-loop hairpin systems consist of two sets of paired short palindromic repeats separated by several unpaired loop sequences, playing a crucial role in physiological processes such as gene expression regulation. Under physiological conditions, RNA can spontaneously fold to form a hairpin stem-loop structure. The ability to reproduce this de novo folding process in simulations is an important indicator of the global rationality of molecular force field parameters. We evaluated the ability of the BSFF1 force field to de novo folding of the r(gcGCAAgc) stem-loop hairpin system using 500 ns copy exchange REMD enhanced sampling simulations. The results are as follows: Figure 4 As shown, the BSFF1 parameters can successfully reproduce all the important pairing interactions on the stem and ring, and the simulated structure is similar to the experimental structure, which further verifies the global rationality of the BSFF1 force field parameters.

[0067] Figure 4 The stem-ring hairpin system was folded from the root, resulting in... Figure 4 It can be seen that BSFF1 can successfully reproduce all the important hydrogen bond pairing interactions on the stem and ring structures with very small deviations from the experimental structures.

[0068] Furthermore, we analyzed the de novo folding path of r(gcGCAAgc) under the BSFF1 parameters, such as... Figure 5 As shown, this result has important guiding value for studying de novo folding of related RNA systems.

[0069] Figure 5 For the stem-ring hairpin system, the folding path is from the head. Figure 5 It is known that the BSFF1 force field can cause RNA to spontaneously fold from its initial disordered structure into a natural stem-loop hairpin structure. This folding pathway has important guiding significance for studying the de novo folding of the corresponding RNA.

[0070] The method for generating base-specific molecular force field parameters (BSFF1) of nucleic acids based on reweighting provided by the embodiment of the present application fully considers the structural property differences of each base and can effectively improve the simulation effect on nucleic acid molecules. The method has reference value for developing more accurate force fields, and the base-specific processing force field development idea is also helpful for the future development of RNA molecular force fields, and is of great help to the research on nucleic acid analysis structure and function, and is conducive to the development of nucleic acid industry such as RNA drugs.

[0071] The basic principles and main features of the present application and the advantages of the present application are shown and described above. Those skilled in the art should understand that the present application is not limited by the above examples, and the above examples and descriptions in the specification are only to illustrate the principles of the present application. Without departing from the spirit and scope of the present application, various changes and improvements can be made to the present application, and these changes and improvements all fall within the scope of the claimed present application. The scope of protection of the present application is defined by the appended claims and their equivalents.

Claims

1. A method for generating base-specific force field parameters for nucleic acid molecules based on a reweighting algorithm, characterized in that, Whether the four-nucleotide conformation generated by molecular dynamics simulation is an incorrect intercalation conformation is used as the basis for specifically adjusting the non-bonding parameters of nucleic acid nucleoside heavy atoms and introducing grid-based energy correction parameters; including the following steps: (1) Initial molecular dynamics simulation: based on nucleic acid force field ff99bsc0χOL3 Parameters for REST2 simulations with solute tempering replica exchange enhancement sampling of the four tetranucleotide systems AAAA, CCCC, GGGG, and UUUU. (2) Trajectory analysis: extract all conformations of the simulation reference copy trajectory, analyze the base stacking order using DSSR software, and write scripts to distinguish whether it is an incorrect intercalation conformation; (3) Re-weighted adjustment of non-bonding parameters: based on whether the conformation is an intercalation structure, the non-bonding parameters of nucleoside heavy atoms on the four bases are specifically adjusted; (4) Second molecular dynamics simulation: based on the adjusted parameters, REST2 simulation with solute tempering copy exchange enhancement sampling is performed on the five four-nucleotide systems of AAAA, CAAU, CCCC, GACC, and UUUU; (5) Trajectory analysis: extract all conformations of the simulation reference copy trajectory, analyze the stacking using DSSR software, and write scripts to distinguish whether it is an incorrect intercalation conformation, and calculate the non-bonding energy between nucleoside heavy atoms; (6) Re-weighted adjustment of CMAP parameters: based on whether the conformation is an intercalation structure, introduce ζ / α main chain dihedral angle grid energy terms through re-weighted algorithm; (7) Obtain a set of base-specific accurate molecular force field parameters and carry out testing.

2. The method of claim 1, wherein the method is based on a reweighting algorithm to generate the force field parameters for the base-specific nucleic acid molecule. Use self-scripted scripts to analyze DSSR software calculation results to determine whether it is an incorrect intercalation conformation.

3. The method for generating base-specific force field parameters for nucleic acid molecules based on a reweighting algorithm according to claim 1 or 2, characterized in that, Using the re-weighted algorithm, based on the proportion of incorrect intercalation conformations, adjust the non-bonding parameters of nucleoside heavy atoms in a base-specific manner.

4. The method of claim 3, wherein the method is based on a reweighting algorithm to generate the force field parameters for the base-specific nucleic acid molecule. Using the re-weighted algorithm, based on the proportion of incorrect intercalation conformations, introduce ζ / α main chain grid energy correction parameters to obtain the final nucleic acid base-specific molecular force field parameters.

Images