A method of developing base-specific force field parameters

Abstract

The application relates to the technical field of cheminformatics, and provides a method for developing base-specific force field parameters, which comprises the following steps: S1, calculating a base pairing energy item; S2, simulating a specific energy item; S3, introducing a molecular force field: adding the pairing energy item BP-fix and BS-CMAP parameters into an existing base-specific force field. By introducing the additional energy items BP-fix and BS-CMAP, the simulation accuracy and sampling efficiency are improved, and the additional energy items introduced have reference value for the development of other RNA force fields; the method is important for the research on the dynamic ensemble structure function mechanism of RNA, the development of RNA targeted drugs, the research and development of RNA therapy, and the like, and is favorable for the development of the nucleic acid medical industry.

Description

Technical Field

[0001] This invention relates to the field of cheminformatics, and particularly to theories applicable to quantum chemistry, molecular mechanics, molecular dynamics, etc. Background Technology

[0002] Ribonucleic acid (RNA) plays a crucial role in many biological processes, including transcription and translation regulation, with RNA structure and dynamics being essential to its functional mechanisms. However, resolving RNA structure and dynamics using traditional experimental methods is extremely difficult. Molecular dynamics simulations (MD), capable of continuous conformational sampling at the atomic scale, have become an important method in RNA research. Current MD simulations still face two major shortcomings: First, although computing power is constantly improving and supercomputing is developing rapidly, current MD simulations are far from sufficient to reproduce many biological processes (sampling efficiency issue); second, the accuracy of existing molecular force fields is insufficient, leading to sampling distortion that fails to reflect the true situation (force field accuracy issue).

[0003] Specifically, the inventors discovered that the RNA molecular force field in the prior art has at least the following problems: First, it seriously underestimates the stability of base pairing; second, it produces a large number of erroneous main chain conformations in the simulation that do not match the experiments. Summary of the Invention

[0004] The purpose of this invention is to provide a method for developing base-specific force field parameters. By introducing the base pairing energy repair term BP-Fix and the base-specific BS-CMAP energy term, a new generation of base-specific nucleic acid molecular force field parameters are generated, thereby improving the accuracy of force field simulation.

[0005] To address the aforementioned technical problems, this invention provides a method for developing base-specific force field parameters, comprising the following steps:

[0006] S1. Calculate the base pairing energy term: Scan the RNA base pairing dimer in solution to calculate the quantum mechanical energy and molecular mechanical energy of the intermolecular distance. Calculate the harmonic confinement energy term based on the quantum mechanical energy and molecular mechanical energy and formulate the pairing energy term BP-fix. RNA base pairing includes AU base pairs, GC base pairs and GU base pairs.

[0007] S2. Simulation calculation of specific energy terms: REST2 simulations with solute tempering were performed on four tetranucleotide systems: AAAA, CCCC, GGGG, and UUUU. Based on the REST2 simulation, the trajectory of the reference copy was extracted, the simulated structure was extracted, and a reweighting method was used to re-weight each conformation based on the convergent conformation ensemble under the existing force field parameters. Then, gamma-delta dihedral angle energy correction was performed on the four tetranucleotide systems respectively, and four sets of base-specific BS-CMAP parameters were generated accordingly.

[0008] S3. Introduce a molecular force field: Add the pairing energy terms BP-fix and BS-CMAP parameters to the existing base-specific force field.

[0009] In step S1, a scan is performed at an interval of 0.05 Å.

[0010] The REST2 simulation uses 16 replicas, with the replica temperature range set from 275 K to 500 K.

[0011] The temperature for different copies is set according to the following formula:

[0012] ,

[0013] In the formula, For the simulated temperature corresponding to copy i, The lowest temperature among all copies. It is the highest temperature among all replicas, and n is the number of replicas.

[0014] The REST2 simulation has 16 replicas with temperatures of 275.00 K, 286.19 K, 297.82 K, 309.93 K, 322.53 K, 335.64 K, 349.29 K, 363.49 K, 378.27 K, 393.65 K, 409.66 K, 426.32 K, 443.65 K, 461.69 K, 480.46 K, and 500.00 K, respectively. The simulation time for each replica is set to 200 nanoseconds, the simulation step size is set to 2 femtoseconds, a frame of result is output every 50 picoseconds, and a replica swap attempt is performed every 1000 steps.

[0015] In step S2, the reweighting calculation is performed using the following formula:

[0016] ,

[0017] In the formula, Indicates average, and They represent the parameters according to the force field. and The ensemble properties of each individual molecular dynamics simulation. and They represent the parameters respectively. and The corresponding potential energy of each conformation is calculated below. The constant is the reciprocal of the product of the Boltzmann constant and the thermodynamic temperature.

[0018] In step S2, the extracted simulated structure is analyzed using DSSR software to determine its base stacking sequence and a self-written script is used to determine whether it is the experimental conformation.

[0019] In step S2, before performing the REST2 simulation, the simulation system is preprocessed using the following steps:

[0020] S2.1 For tetranucleotides in the standard conformation, the pdb2gmx command of Gromacs software is used to perform topological file conversion, and the Tip3p water model and resistance ions are introduced to simulate physiological conditions.

[0021] S2.2 Use the steepest descent method to perform energy minimization optimization on the above system in a maximum of 4000 steps;

[0022] S2.3. The above system is heated to 275K by 5000 NVT ensemble heating and 25000 NPT ensemble heating.

[0023] In step S1, the distance between the atomic pairs is defined as follows: AU base pair is A@N1-U@H3, GC base pair is G@H1-C@N3, and GU base pair is G@H1-U@O2.

[0024] In step S1, quantum mechanical energy and molecular mechanical energy are calculated simultaneously under the condition of implicit solvent. The quantum mechanical energy is calculated using Gaussian09 at the M05-2X / 6-311G** / SMD theoretical level, and the molecular mechanical energy is calculated using OL3 and BSFF1 force fields in Amber, using the GBneck2 implicit solvent model and the mbondi3 radius set for polar solvation and SASA-based nonpolar solvation.

[0025] Compared with existing technologies, this invention improves the accuracy and sampling efficiency of simulation by introducing additional energy terms BP-Fix and BS-CMAP. At the same time, the introduced additional energy terms are also of reference value for the development of other RNA force fields. It is of great help to study the dynamic ensemble structure and function mechanism of RNA, develop RNA-targeted drugs, and develop RNA therapies, which is conducive to the development of the nucleic acid medical industry.

[0026] The above description is merely an overview of the technical solution of the present invention. In order to better understand the technical means of the present invention and to implement it in accordance with the contents of the specification, and in order to make the above and other objects, features and advantages of the present invention more apparent and understandable, specific embodiments of the present invention are described below. Attached Figure Description

[0027] One or more embodiments are illustrated by way of example with the corresponding pictures in the accompanying drawings. These illustrations do not constitute a limitation on the embodiments. Elements with the same reference numerals in the drawings are denoted as similar elements. Unless otherwise stated, the figures in the drawings do not constitute a limitation on scale.

[0028] Figure 1 This is a flowchart illustrating at least one embodiment of the present invention;

[0029] Figure 2 This is a graph showing the test results of a tetranucleotide system according to an embodiment of the present invention;

[0030] Figure 3 This is a diagram showing the result of folding the four-ring hair clip system from the beginning according to an embodiment of the present invention. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of this invention clearer, the various embodiments of this invention will be described in detail below with reference to the accompanying drawings. However, those skilled in the art will understand that many technical details have been provided in the various embodiments of this invention to facilitate a better understanding of this application. However, the technical solutions claimed in this application can be implemented even without these technical details and with various variations and modifications based on the following embodiments. The division of the following embodiments is for ease of description and should not constitute any limitation on the specific implementation of this invention. The embodiments can be combined with and referenced by each other without contradiction.

[0032] The first embodiment of the present invention relates to a method for developing base-specific force field parameters, the process of which is as follows: Figure 1 As shown, it includes the following steps:

[0033] S1. Calculate the base pairing energy term: Scan the RNA base pairing dimer in solution to calculate the quantum mechanical energy and molecular mechanical energy of the intermolecular distance. Calculate the harmonic confinement energy term based on the quantum mechanical energy and molecular mechanical energy and formulate the pairing energy term BP-fix. RNA base pairing includes AU base pairs, GC base pairs and GU base pairs.

[0034] S2. Simulation calculation of specific energy terms: REST2 simulations with solute tempering were performed on four tetranucleotide systems: AAAA, CCCC, GGGG, and UUUU. Based on the REST2 simulation, the trajectory of the reference copy was extracted, the simulated structure was extracted, and a reweighting method was used to re-weight each conformation based on the convergent conformation ensemble under the existing force field parameters. Then, gamma-delta dihedral angle energy correction was performed on the four tetranucleotide systems respectively, and four sets of base-specific BS-CMAP parameters were generated accordingly.

[0035] S3. Introduce a molecular force field: Add the pairing energy terms BP-fix and BS-CMAP parameters to the existing base-specific force field.

[0036] Therefore, this invention, on the one hand, uses quantum mechanics and molecular mechanics to quantify the problem of excessively high base pairing energy and low stability in existing molecular force fields, and rationally designs by using the difference between the two, introducing the base pairing energy repair term BP-Fix; on the other hand, it uses simulation calculation analysis on a tetranucleotide system, combined with a reweighted algorithm, to introduce a base-specific gamma-delta main chain dihedral energy term BS-CMAP based on whether the simulated structure is the experimental conformation, ultimately achieving the goal of a new generation of base-specific force field parameters.

[0037] The steps of the various methods described above are only for clarity. In practice, they can be combined into one step or some steps can be split into multiple steps. As long as they include the same logical relationship, they are all within the scope of protection of this patent. Adding insignificant modifications or introducing insignificant designs to the algorithm or process, but without changing the core design of the algorithm and process, are also within the scope of protection of this patent.

[0038] The second embodiment of the present invention is largely the same as the first embodiment, except that step S1 is further optimized, in which scanning is performed at an interval of 0.05 Å.

[0039] Furthermore, in step S1, the distance between the atomic pairs is defined as follows: AU base pair is A@N1-U@H3, GC base pair is G@H1-C@N3, and GU base pair is G@H1-U@O2.

[0040] Furthermore, in step S1, the quantum mechanical energy and molecular mechanical energy are calculated simultaneously under the condition of implicit solvent. The quantum mechanical energy is calculated using Gaussian09 at the M05-2X / 6-311G** / SMD theoretical level, and the molecular mechanical energy is calculated using OL3 and BSFF1 force fields in Amber, using the GBneck2 implicit solvent model and the mbondi3 radius set for polar solvation and SASA-based nonpolar solvation.

[0041] The third embodiment of the present invention is largely the same as the first embodiment, except that step S2 is further optimized. The REST2 simulation uses 16 replicas, and the temperature range of the replicas is set to 275 K to 500 K.

[0042] Furthermore, the temperature for different replicas is set according to the following formula:

[0043] ,

[0044] In the formula, For the simulated temperature corresponding to copy i, The lowest temperature among all copies. It is the highest temperature among all replicas, and n is the number of replicas.

[0045] Preferably, the REST2 simulation has 16 replicas with temperatures of 275.00 K, 286.19 K, 297.82 K, 309.93 K, 322.53 K, 335.64 K, 349.29 K, 363.49 K, 378.27 K, 393.65 K, 409.66 K, 426.32 K, 443.65 K, 461.69 K, 480.46 K, and 500.00 K, respectively. The simulation time for each replica is set to 200 nanoseconds, the simulation step size is set to 2 femtoseconds, one frame of result is output every 50 picoseconds, and a replica swap attempt is performed every 1000 steps.

[0046] Furthermore, in step S2, the reweighting calculation is performed using the following formula:

[0047] ,

[0048] In the formula, Indicates average, and They represent the force field parameters respectively. and The ensemble properties of each individual molecular dynamics simulation. and They represent the parameters respectively. and The corresponding potential energy of each conformation is calculated below. The constant is the reciprocal of the product of the Boltzmann constant and the thermodynamic temperature.

[0049] Furthermore, in step S2, the extracted simulated structure is analyzed using DSSR software to determine its base stacking sequence and a self-written script is used to determine whether it is the experimental conformation.

[0050] Furthermore, in step S2, before performing the REST2 simulation, the simulation system is preprocessed using the following steps:

[0051] S2.1 For tetranucleotides in the standard conformation, the pdb2gmx command of Gromacs software is used to perform topological file conversion, and the Tip3p water model and resistance ions are introduced to simulate physiological conditions.

[0052] S2.2 Use the steepest descent method to perform energy minimization optimization on the above system in a maximum of 4000 steps;

[0053] S2.3. The above system is heated to 275K by 5000 NVT ensemble heating and 25000 NPT ensemble heating.

[0054] The fourth embodiment of the present invention, in conjunction with the above embodiments, adopts the following steps:

[0055] (1) Quantum mechanics and molecular mechanics are used to calculate the base pairing energy term and introduce an additional energy term BP-Fix.

[0056] To quantitatively assess the energy bias of the force field on base pairing stability, intermolecular distance (QM) and (MM) energies of RNA base pair dimers were calculated in solution. In addition to the two Watson-Crick base pairs (AU and GC), the most common non-Watson-Crick base pair in RNA, the GU swing base pair (GU base pair), was also introduced. Due to its structural stability and important biological function, it generates base pair dimer conformations at different distances. The distances were defined as follows: AU base pairs were A@N1-U@H3, GC base pairs were G@H1-C@N3, and GU base pairs were G@H1-U@O2. Intermolecular distances between base pairs ranged from 1 Å to 6 Å, with scan intervals of 0.05 Å. Since the BP-fix term should include any errors in base pair interactions with aqueous solvents in the MD simulation, solvent effects needed to be considered in the calculations; therefore, both QM and MM energies were calculated simultaneously under conditions containing implicit solvents. This strategy has been successfully used in RNA and protein force field parameterization to address the inconsistency between the gas and solution phases. QM calculations were performed using Gaussian09 at the M05-2X / 6-311G** / SMD theoretical level. MM calculations were performed in Amber using the OL3 and BSFF1 force fields, employing the GBneck2 implicit solvent model (igb=8) and the mbondi3 radius set for polar solvation and SASA-based nonpolar solvation. Based on the energy difference between the harmonic quantum mechanical (QM) calculations and the molecular mechanical (MM) force fields, the harmonic confinement energy terms were precisely configured to formulate the final BP-fix term.

[0057] (2) Initial molecular dynamics simulation of tetranucleotides.

[0058] To improve sampling accuracy and ensure sufficient simulation, an enhanced sampling method with solute tempering and copy exchange (REST2) was first used. This method aimed to converge the simulation trajectory as much as possible using limited computational resources. REST2 simulations were performed on four tetranucleotide systems: AAAA, CCCC, GGGG, and UUUU using GROMACS 2018.8 and plumed 2.6.2 software. Preprocessing of the simulation systems included: first, topology file conversion for tetranucleotides with standard conformations was performed using the pdb2gmx command in Gromacs software, introducing the Tip3p water model and resistance ions to simulate physiological conditions; second, the steepest descent method was used to perform energy minimization optimization on the above systems for up to 4000 steps; and third, the above systems were subjected to 5000 steps of NVT ensemble heating and 25000 steps of NPT ensemble heating to 275K. Throughout the process, all covalent bonds connected to hydrogen atoms were constrained using the LINCS algorithm, while long-range electrostatic interactions were calculated using the Particle Mesh Ewald (PME) algorithm. The cutoff values ​​for Lennard-Jones interactions and electrostatic interactions were both set to 10 Å (1 nm). After the above treatment, a REST2 simulation with solute tempering was performed. The simulation used 16 exchanges, and the replica temperatures were set to a range of 275 K to 500 K. Based on the principle of REST2, the temperatures for different replicas were set according to formula (1):

[0059] ……………………………(1)

[0060] in, This represents the simulated temperature of the corresponding copy i. It is the lowest temperature among all replicas, that is, the temperature of reference replica 0. is the highest temperature among all replicas, and n is the number of replicas. According to formula (1), we obtain the temperatures of the 16 replicas as 275.00 K, 286.19 K, 297.82 K, 309.93 K, 322.53 K, 335.64 K, 349.29 K, 363.49 K, 378.27 K, 393.65 K, 409.66 K, 426.32 K, 443.65 K, 461.69 K, 480.46 K, and 500.00 K, respectively. The simulation time for each replica is set to 200 nanoseconds, the simulation step size is set to 2 femtoseconds, one frame of result is output every 50 picoseconds, and a replica swap attempt is performed every 1000 steps (2 picoseconds).

[0061] (3) Trajectory analysis and feature calculation.

[0062] For the initial molecular dynamics simulation trajectory, the trajectory of the reference copy (the copy at the lowest temperature) was extracted, the simulated structure was extracted, and its base packing sequence was analyzed using DSSR software. A self-written script was used to determine whether it was the experimental conformation. At the same time, its dihedral angle information was calculated and saved.

[0063] (4) Introduce a base-specific BS-CMAP energy term through a reweighting algorithm.

[0064] The reweighting method can reweight each conformation under new force field parameters based on the convergent conformation ensemble under the existing force field parameters. Thus, the new ensemble and the proportion of each conformation can be obtained directly without resimulating, and then the average property can be obtained, i.e., formula (2).

[0065] …………………(2)

[0066] in, Indicates average, and They represent the parameters according to the force field. and The ensemble properties of each individual molecular dynamics simulation. and They represent the parameters respectively and The corresponding potential energy of each conformation is calculated below. The constant is the reciprocal of the product of the Boltzmann constant and the thermodynamic temperature.

[0067] The reweighting method can quantitatively link force field parameters with simulated properties, thereby obtaining simulated properties under new parameters through simulation with existing parameters, and thus evaluating the rationality of the new force field parameters. In this invention, we use the reweighting method to perform gamma-delta dihedral energy correction on four tetranucleotide systems based on whether each frame structure is the experimental conformation and its corresponding dihedral angle information, thereby generating four sets of base-specific BS-CMAP parameters.

[0068] (5) A self-written script introduces the above-mentioned BP-Fix and BS-CMAP additional energy terms into the molecular force field.

[0069] The additional energy terms of BP-Fix and BS-CMAP obtained in steps 1 and 4 are added to the existing molecular force field through a script, thus finally obtaining the next generation of base-specific force field parameters (BSFF2).

[0070] (6) Molecular dynamics simulation test.

[0071] The REST2 simulation method mentioned in step 1 was used to test the simulation effect of BSF2 parameters in the tetranucleotide system and compare it with the ff99bsc0χOL3 force field. The simulation time and conditions were kept consistent with those in step 1. Conventional molecular dynamics simulation method was used to simulate de novo folding of the hairpin stem-loop system using the BSFF2 force field. The temperature was set to 308K and the simulation time was 1000 nanoseconds.

[0072] (7) Post-molecular dynamics simulation analysis.

[0073] After simulation, the trajectory was exported as a pdb structure using the trjconv module of Gromacs, and the root mean square deviation (RMSD) was calculated using the cpptraj module of AMBER and the DBSCAN module was used for structural clustering.

[0074] Example 1

[0075] Based on the fourth embodiment described above, the accuracy of the BSFF2 force field is evaluated by performing enhanced molecular dynamics simulations on tetranucleotides. At the same time, the BSFF2 force field can fold the hairpin stem-loop system from scratch in a short time through conventional molecular dynamics simulations, thereby further verifying the global rationality of the parameters and the sampling efficiency.

[0076] 1. Tetranucleotide:

[0077] The tetranucleotide system is small and has abundant experimental values, making it the gold standard system for testing the performance of RNA molecular force fields. Previous RNA force fields were prone to producing erroneous intercalation conformations when simulating tetranucleotides, resulting in low simulation accuracy. In this invention, we combined the TIP3P water model with the most commonly used nucleic acid force field ff99bsc0χOL3 and our developed BSFF2 force field to perform REST2 enhanced sampling on five tetranucleotide systems: AAAA, CCCC, UUUU, CAAU, and GACC, and conducted structural cluster analysis. The results are as follows: Figure 2 As shown, it can be found that the ff99bsc0χOL3 force field will produce a large number of erroneous intercalation structures. The BSFF2 simulation significantly improved the simulation sampling, not only successfully reducing the occurrence of erroneous intercalation conformations, but also increasing the frequency of reasonable experimental conformations, thus greatly improving the simulation accuracy.

[0078] 2. The stem ring hair clip system is folded from the top:

[0079] RNA stem-loop hairpin systems, especially tetraloop systems, are important structural modules of RNA, playing crucial roles in physiological functions such as protein expression and gene regulation. In vivo, RNA stem-loop hairpins can spontaneously fold, a highly complex process. Reproducing this de novo folding process in simulations demands high global accuracy and sampling efficiency of the molecular force field. For a tetraloop structure with a PDB ID of 8CLR, we performed a 1-microsecond conventional molecular dynamics simulation using a BSFF2 force field at the experimental temperature, successfully achieving de novo folding. Figure 3 As shown, the structure exhibits an extended configuration during the first 25 ns of simulation, and partially folds and undergoes terminal base pairing at 50 ns. The initial folding simulation was achieved within the first 100 ns. During the remaining 900 ns of simulation, BSFF2 fully sampled various conformations of the tetracyclic ring, with the closest RMSD to the experimental structure being 0.7 Å. This further validates the overall accuracy of BSFF2, and compared to previous force fields, BSFF2 demonstrates extremely high simulation sampling efficiency.

[0080] Those skilled in the art will understand that the above embodiments can be modified in form and detail in practical applications without departing from the spirit and scope of the invention.

Claims

1. A method for developing base-specific force field parameters, characterized in that, Includes the following steps: S1. Calculate the base pairing energy term: Scan the RNA base pairing dimer in solution to calculate the quantum mechanical energy and molecular mechanical energy of the intermolecular distance. Calculate the harmonic confinement energy term based on the quantum mechanical energy and molecular mechanical energy and formulate the pairing energy term BP-fix. RNA base pairing includes AU base pairs, GC base pairs and GU base pairs. S2. Simulation calculation of specific energy terms: REST2 simulations with solute tempering were performed on four tetranucleotide systems: AAAA, CCCC, GGGG, and UUUU. Based on the REST2 simulation, the trajectory of the reference copy was extracted, the simulated structure was extracted, and a reweighting method was used to re-weight each conformation based on the convergent conformation ensemble under the existing force field parameters. Then, gamma-delta dihedral angle energy correction was performed on the four tetranucleotide systems respectively, and four sets of base-specific BS-CMAP parameters were generated accordingly. S3. Introduce a molecular force field: Add the pairing energy terms BP-fix and BS-CMAP parameters to the existing base-specific force field.

2. The method for developing base-specific force field parameters according to claim 1, characterized in that, In step S1, a scan is performed at an interval of 0.05 Å.

3. The method for developing base-specific force field parameters according to claim 1, characterized in that, The REST2 simulation uses 16 replicas, with the replica temperature range set from 275 K to 500 K.

4. The method for developing base-specific force field parameters according to claim 3, characterized in that, The temperature for different copies is set according to the following formula: ， In the formula, For the simulated temperature corresponding to copy i, The lowest temperature among all copies. It is the highest temperature among all replicas, and n is the number of replicas.

5. The method for developing base-specific force field parameters according to any one of claims 1 or 3, characterized in that, The REST2 simulation has 16 replicas with temperatures of 275.00 K, 286.19 K, 297.82 K, 309.93 K, 322.53 K, 335.64 K, 349.29 K, 363.49 K, 378.27 K, 393.65 K, 409.66 K, 426.32 K, 443.65 K, 461.69 K, 480.46 K, and 500.00 K, respectively. The simulation time for each replica is set to 200 nanoseconds, the simulation step size is set to 2 femtoseconds, a frame of result is output every 50 picoseconds, and a replica swap attempt is performed every 1000 steps.

6. The method for developing base-specific force field parameters according to claim 1, characterized in that, In step S2, the reweighting calculation is performed using the following formula: ， In the formula, Indicates average, and They represent the parameters according to the force field. and The ensemble properties of each individual molecular dynamics simulation. and They represent the parameters respectively. and The corresponding potential energy of each conformation is calculated below. The constant is the reciprocal of the product of the Boltzmann constant and the thermodynamic temperature.

7. The method for developing base-specific force field parameters according to claim 1, characterized in that, In step S2, the extracted simulated structure is analyzed using DSSR software to determine its base stacking sequence and a self-written script is used to determine whether it is the experimental conformation.

8. The method for developing base-specific force field parameters according to claim 1, characterized in that, In step S1, the distance between the atomic pairs is defined as follows: AU base pair is A@N1-U@H3, GC base pair is G@H1-C@N3, and GU base pair is G@H1-U@O2.

9. The method for developing base-specific force field parameters according to claim 1, characterized in that, In step S1, quantum mechanical energy and molecular mechanical energy are calculated simultaneously under the condition of implicit solvent. The quantum mechanical energy is calculated using Gaussian09 at the M05-2X / 6-311G** / SMD theoretical level, and the molecular mechanical energy is calculated using OL3 and BSFF1 force fields in Amber, using the GBneck2 implicit solvent model and the mbondi3 radius set for polar solvation and SASA-based nonpolar solvation.

Images