Extensive Evaluation of Force Fields for G-Quadruplexes
G-Quadruplexes (GQs), folded by guanine-rich nucleic acid sequences, involve in gene expression processes and closely associated with the formation of tumors. So far, GQ has drawn widespread attention for its notable application of serving as potential anti-cancer target. Recently, theoretical studies for GQs have achieved significant progress, most of which are inseparable from molecular dynamics (MD) simulation. As a necessary tool to explore dynamics behavior of molecules, MD simulations strictly depend on force field parameters, which is a sticking point to obtain accurate results. Currently, many force fields for nucleic acids have been developed, but none of them have been accepted widely for GQs. In this paper, we selected five popular force fields, which are parmbsc0, parmbsc1, OL15, Drude2017 and AMOEBANUC17, and conducted explicit-solvent MD simulations on two DNA GQs respectively. We evaluated these force fields from many aspects in detail. Meanwhile, we compared conformational energy using quantum chemistry calculations.With the comprehensive evaluation, Drude2017 achieved better description for GQs, which we suggest that using Drude2017 force field should be taken into account first when investigating GQs by MD simulation.
INTRODUCTION<!>Starting systems<!>Simulation protocol<!>Conformational energy calculation<!>Stability of the GQs<!>Distance between channel ions<!>Stability of Hoogsten hydrogen bonds<!>Compare dihedrals with experiment<!>Conformational energy<!>CONCLUSION AND PERSPECTIVE
<p>G-Quadruplexes (GQs) 1 are folded by guanine-rich nucleic acid sequences. Every four guanine constitutes a planar square called tetrad and interact with neighbor guanines by hoogsten hydrogen bonds. The curved nucleic acid strands prompt each guanine tetrad to stack layer by layer with a linear cavity formed. Metal cations coordinate with carbonyl oxygens in the cavity to neutralize negative charges of the system for structural stability. [2][3][4] GQs exist in chromosomal telomeres and promoter region, which play a vital important role in regulating gene expression like replication, transcription, translation, and thus fundamentally affect cellular proliferation especially in cancer cells. [5][6][7] In recent years, there have been increasing appeal of GQ-targeting therapeutics, for this is a promising pathway to treat cancer and other human diseases. [8][9][10][11][12][13] In addition, GQs are widely used to construct biosensors, catalyze reactions and design novel materials owing to its unique properties of ionic affinity and easy modification. [14][15] Thus, it is necessary to understand the structure and function of GQ thoroughly. However, presently we still have a limited knowledge to clarify the correlation between GQ conformation and function exactly, which motivates intense exploration of GQs.</p><p>Since the GQ was first identified in the 1960s, 16 considerable endeavors have been made to study its formation, physiological mechanism and application. 17 High-resolution characterization techniques such as nuclear magnetic resonance, X-ray crystallography, and cryogenic electron microscopy are capable of investigating GQ structures. Meanwhile, theoretical computation methods based on quantum mechanics (QM) and molecular mechanics (MM) provide deep insight into kinetic and thermodynamic properties. Given that some minute changes, for example, sequences, channel ions, base orientations and folding modes contribute to structural polymorphism of GQ, [18][19] theoretical study especially Molecular Dynamics (MD) simulation [20][21] is increasingly an indispensable tool to study GQs. MD simulation displays how molecules respond to controllable variables by modeling them in real physiological environment at full atomic and femtosecond level, which is an effective way to observe molecular motion.</p><p>Noteworthily, reliable results of MD simulation highly depend on accurate force fields.</p><p>Currently there are different kind of force fields designed for canonical DNA or RNA. Over the past few years, some of them have been actually applied to conduct MD simulation of GQs, which achieved great success in interpreting GQs dynamics behavior, including their folding pathways. [22][23][24][25] However, several flaws existing in the simulation have been raised subsequently, for example, bifurcated hydrogen bonds 22 observed in simulations but not discovered in the crystal structures.</p><p>Such inconsistence was attributed to intrinsic approximation in pair-additive force fields. 26 Meanwhile, many reviews revealed the overestimation of short-range interaction in the MD simulation. [27][28][29] The pair-additive force fields adopt fixed atomic partial charges so that induced polarization effect is not considered. This defect can be magnified as GQs are flexible and highly electronegative in real environment. Therefore, the pair-additive force fields cause an average effect of charges and the corresponding results should be discussed cautiously.</p><p>For years, there have been a remarkable improvement of force fields to describe GQs, particularly polarizable force field applied in MD simulation. At first, force fields parm94 30 and parm98/99 31 were used broadly to simulate GQs. Whereas in 2007, Pérez et al. refined parm99 into parmbsc0 force field by modifying α/γ torsional term, 32 subsequently followed by OL15 33 and parmbsc1 34 refined force fields. They are more recently dominant Amber force fields to implement MD simulation for DNA, also for GQ. [35][36][37] Moreover, CHARMM force fields are likewise employed extensively, yet of which CHARMM27 has been validated to be far from satisfactory in MD simulation of GQ. 27 In 2014, Savelyev and MacKerell released a CHARMM polarizable force field termed Drude2013 for DNA based on the classical Drude oscillator model. 38 Its recent modification Drude2017 was claimed that it outperforms CHARMM36 force field in the simulation of GQ. [39][40][41][42] Besides, with induced polarization effect included, Song et al. developed polarized nucleic acidspecific charge, namely PNC force field, and obtained a significantly improved description of GQ compared with Amber ff10 force field. 26 The force field, atomic multipole optimized energetics for biomolecular applications (AMOEBA) also incorporates polarization effect. 43 It is capable of describing nucleic acids, but so far it has not been applied for GQ in relevant publications.</p><p>Furthermore, some studies also accessed the influence of water model and ion parameters, which facilitates a more reliable MD simulation for GQ. [44][45] As there are many force fields available to choose from to simulate GQs, but none of them are widely accepted. Hence, in this work, we conducted a systematic test and evaluation of force fields for GQ. We selected five force fields that are parmbsc0, parmbsc1, OL15, Drude2017 and AMOEBANUC17, for they are most recently used to study GQ. These force fields were respectively applied to simulate two DNA GQs containing K + and Na + in channel. Since the two kinds of ions are physiology-related and possess different stability of GQ, it is appropriate to utilize them to test the performance of force fields respectively. To make a persuasive comparison among these force fields, we calculated conformational energy of modified GQs using MM and QM. Our study would highlight the selection of force fields to simulate GQ.</p><!><p>The starting potassium-containing GQ was gained from protein data bank (PDB ID 1JRN). It is a crystal structure of Oxytricha nova telomeric DNA GQ with the sequence of d[G4T4G4]2, 46 including four chains formed into two separate and identical GQs. We removed one of them as we only need one GQ. The rest GQ holds five potassium ions located in channel with four tetrads and two diagonal loops. (Fig. 1A) Guanine bases in the same strand adopt alternating syn-anti glycosidic bond conformations. The starting sodium-containing GQ was also obtained from protein data bank (PDB ID 1JB7). 47 It is a part of Oxytricha nova telomeric protein-DNA complex characterized by x-ray diffraction. Since GQ is the system of interest and there is a certain distance between the protein and DNA GQ, we removed protein from the complex. Four sodium ions are situated in plane of four tetrads. (Fig. 1B) The two GQs conformations are completely identical except for ion species. Also, we removed the water in initial crystal structures for both GQs.</p><!><p>As the process of MD simulation for both GQs is almost the same, here we first take potassiumcontaining GQ for example. All the MD simulations were carried out in explicit solvent with periodic boundary conditions on GPU. "Drude", "AMOEBA", "bsc0" and "bsc1" are the abbreviation for Drude2017, AMOEBANUC17, parmbsc0 and parmbsc1 force field.</p><p>For Amber force fields, we used AMBER18 software package, 48 of which tLEaP program was used to generate topologies and coordinates of system under OL15, bsc1, and bsc0 force field successively. For potassium-containing GQ, we solvated it in a truncated octahedral box of SPC/E water with the buffer of 12Å. The water model was suggested by a recent study. 44 To neutralize negative charge of the system, 17 potassium ions that were Joung and Cheatham parameters were added into solution as counterions. 49 The prepared system was first minimized for 50000 steps with 200kcal/(mol•Å 2 ) restraints on GQ and five channel potassium ions, which only allowed water molecules and counterions to move freely. This consisted of 25000-step steepest descent followed by 25000-step conjugated gradient. Next, the system was heated gradually from 0K to 300K for 100ps with GQ and channel potassium ions restrained at 200kcal/(mol•Å 2 ). Then an unrestrained equilibration was performed for 100ps in NPT ensemble at constant 1 atm and 300K. Finally, a 200ns production was executed with the same condition as equilibration. The temperature was regulated via Langevin thermostat. Particle Mesh Ewald (PME) was applied to treat long-range electrostatic interactions. 50 A cutoff of 10Å accounted for non-bonded interactions. Hydrogen atoms were constrained by SHAKE algorithm 51 and the time step is 2 fs. Cpptraj 52 in Ambertools was used to analyze trajectories.</p><p>We used CHARMM-GUI 53 to construct solvated potassium-containing GQ under Drude2017 force field. The GQ was placed in an octahedral water box filled with 0.15M KCl via Solution Builder module in CHARMM-GUI. Then topology and coordinate files of the solvated GQ were put into Drude Prepper to add Drude particles and lone pairs of non-hydrogen atoms. The final system was used to conducted MD simulation with OpenMM. 54 The system was minimized followed by NPT equilibration at 1 atm and 303.15 K through Langevin dynamics. The temperature of Drude particles was set to 1 K. Monte Carlo barostat regulated the pressure of system. The integration time step was 1 fs. PME was used to calculate electrostatic interactions and van der Waals forces were switched to zero from 10Å to 12Å. Finally, a 200ns production was performed, and parameters were set the same as equilibration.</p><p>For AMOEBA force field, we used Tinker-OpenMM software 55 to carry out MD simulation.</p><p>The force field applied exactly was AMOEBANUC17 in Tinker8.7 version. 56 To construct a water box, we first built a water molecule using Avogadro software. 57 The water molecule was put into Tinker and a periodic rectangular box was created with the size of 52Å×52Å×52Å through xyzedit program in Tinker. The potassium-containing GQ was soaked in the center of box and 17 potassium ions were placed in box to ensure electric neutrality of the system. The prepared system was minimized using steepest descent gradient optimization at the final RMS gradient of 0.1. During minimization, 50 kcal/(mol•Å 2 ) position restraint was exerted on the GQ and bound ions. Then the same position restraint used in heating process for 500ps in NVT ensemble. The heated system was performed equilibration for 1ns in NPT ensemble at 300K and 1 atm without restraints. In simulation, Andersen thermostat was switched on and a velocity Verlet integrator was adopted. The time step was 1 fs. PME was applied with space cutoff of 7.0Å. The van der Waals cutoff was 12.0Å. The pressure was controlled by Monte Carlo barostat. Likewise, the final production was performed for 200ns.</p><p>For sodium-containing GQ, all the treatment in MD simulation was the same as potassiumcontaining GQ, except for two points. One is that the counterions in solvated sodium-containing GQ were sodium ions, the other is that a rectangular water box was created when using Drude force field.</p><!><p>To compare the performance of the five force fields persuasively, we calculated conformational energy of modified GQs based on MM and QM respectively. From obtained potassium-containing GQ trajectories under Drude2017 force field, we extracted 20 frames that represented 20 different GQ structures (solvent and Drude particles were discarded). Two loops, tetrad 3, tetrad 4, and two bottom potassium ions of each GQ were removed. To ensure the integrity of each GQ structure, corresponding hydrogen atoms were added to terminal chains. Each GQ has two layers of guanine quartets and three potassium ions. (Fig. 2) For MM calculation, each GQ was implemented a very short MD simulation (10 -3 fs) in vacuum with above mentioned Amber force fields respectively at 1K. The conformational energy can be obtained in output files. For AMOEBA force field, it can be directly given using Tinker. For Drude force field, we adopted CHARMM package (version 45b1) 58 to calculate conformational energy of these GQs. We calculated single point energy of each GQ without any optimization using ORCA software (Version 4.1.2). 59 We used B3LYP functional with def2-TZVP(-f) [60][61] basis set and London dispersion correction (D3) 62 with Becke−Johnson (BJ) damping. 63 Resolution of the identity (RI) approximation was used in calculation.</p><p>Fig. 2 The modified potassium-containing GQ. Potassium ions are colored in purple.</p><!><p>Every single MD simulation was up to 200ns. In our work, it is enough time to compare the performance of each force field in the simulations. We first examined the stability of GQs indicated by root mean square deviations (RMSDs) of backbone atoms (C3', C4', C5', O3', O5', P). As shown in Fig. 3, the RMSD are all about within 2.0 Å for potassium-containing GQ in five force fields. In Drude and OL15 force fields, the RMSDs are almost have identical distribution, and they have the lowest value. While for sodium-containing GQ, it has slightly higher RMSD in two polarizable force fields and presents close RMSD in three Amber force fields. Both GQs in AMOEBA force field have the highest RMSD, especially two peaks of RMSD frequency in sodium-containing GQ.</p><p>It implies that both GQs are more flexible in AMOEBA force field. However, the GQs in these force fields are conformationally stable as RMSDs converge at relatively small values, which is qualified to conduct further analysis in detail.</p><!><p>The cations located in GQ channel determine the property of the whole structure to a large extent, particularly its stability. Many studies have revealed stability ranking of conformation for different cations in GQ channel. [64][65][66] Because potassium ion and sodium ion have different size of radius, they have own preferred positions in GQ stem. 41,67 Potassium ions are mainly located between two tetrads, while sodium ions are in the plane of tetrads, just like the two GQs (1JRN and 1JB7) used in our work. (Fig. 1) However, cations in two ends of channel especially in loop region, extremely easily escape from GQ stem in diverse classical MD simulations. This phenomenon has been ascribed to the natural deficiency of force fields proved by QM calculation. [27][28] Thus, in our work we monitored the distance of two adjacent ions at the upper and lower ends of GQ channel respectively. As displayed in Fig. 4, it is noteworthy that none of channel ions escape from channel for both GQs throughout the simulation using Drude force field, which shows the finest and most stable distance of ions. However, for simulations under other four force fields in potassium-containing GQ, at least one end of channel ions went into the solution within 1 ns and never went back to GQ stem. Consequently, there were at most three ions left in the channel. Despite that K27 did not went out of GQ stem in bsc1 force field, the distance of K25-27 has been far beyond normal range. 67 Nevertheless, there is some difference in sodium-containing GQ. All sodium ions within GQ stem were retained in MD simulations using Amber force fields, which the frequency of distance almost overlaps. This is because that pair-addictive force fields have difficulty describing the behavior of channel ions in loop region, whereas no ions exist in the loops of sodium-containing GQ. 27 In AMOEBA force field, both ends of channel ions were lost for potassium-containing GQ.</p><p>Although only one end of channel ion escaped from sodium-containing GQ, the distance of Na25-27 was large. The Na25 appeared very active and even ran out of loop region at the farthest distance of Na25-27 up to 18.27Å, but next moment it went back to approach Na27 rapidly. Obviously, the interaction between sodium ions and GQ stem is not captured correctly. By contrast, Drude force field is of strongly capable of characterizing ion-solute interaction for both GQs.</p><!><p>Adjacent guanine in a tetrad forms two types of circular hoogsten hydrogen bond, which are N1-O6 and N2-N7 respectively. Hoogsten hydrogen bonds play a critical role in maintaining GQ conformational stability. Each tetrad forms four N1-O6 and four N2-N7 hydrogen bonds. We inspected the average distance of the two types of hydrogen bonds as presented in Fig. 5 and Fig. 6.</p><p>For a solid hydrogen bond, its distance should not exceed 3.5 Å at most. In this aspect, Drude force field is outstanding. For both GQs, there were shortest hoogsten hydrogen bonds in Drude force field, which implied stable and strong hydrogen bond interaction. In other force fields, the distance of some N1-O6 and N2-N7 have exceeded 4Å, which obviously cannot form reliable hydrogen bonds. Generally, there is narrow difference in the performance of Amber and AMOEBA force fields in terms of the distance of hoogsten hydrogen bonds. Moreover, for sodium-containing GQ, the frequency of distance distribution almost overlaps completely.  To evaluate the formation of hydrogen bonds in MD simulation in detail, we compared the occupation ratio of each hydrogen bond during MD simulation. (Fig. 7) Likewise, no mater which GQ and which type of hoogsten hydrogen bond, the occupation ratio is always close to percent 100 in Drude force field, much higher than other force fields. It demonstrates longer lifetime of hydrogen bonds in MD simulation. Meanwhile, the angles of hydrogen bonds present more concentrated distributions using Drude force field for both GQs. (Fig. 8) However, the performance of Amber and AMOEBA force fields is similar, which is knotty to make a comparison. Since one tetrad forms four N1-O6 and four N2-N7 hydrogen bonds respectively, 32 hydrogen bonds in total should exist in both GQs if force field is excellent enough. Thus, we examined the number of hydrogen bonds in each frame and presented them with quantity distribution. (Fig. 9) As excepted, the number of hydrogen bonds varies from 25 to 32 and most of them are close to 32 in Drude force field for both GQs. However, the quantity has a wider distribution, mainly from 15 to 20 in Amber and AMOEBA force fields.   In addition, the bifurcated hydrogen bond that is N1-N7 has been observed in previous MD simulations. 22 This kind of hydrogen bond is regarded as an artifact due to the inaccuracy of force field. Therefore, we inspected the existence of bifurcated hydrogen bond in our MD simulations.</p><p>(Table 1) For Drude force field, encouragingly, we found that none of such hydrogen bonds formed throughout the MD simulation in potassium-containing GQ, whereas they existed momentarily in sodium-containing GQ. For Amber force fields, the occupation ratios of bifurcated hydrogen bonds are more than percent 20 in sodium-containing GQ. For AMOEBA force field, bifurcated hydrogen bonds are formed with a short lifetime in both GQs. By comparison, Drude still outperforms other force fields.</p><!><p>To further examine the performance of these five force fields, comparison with experiment is needed.</p><p>Both GQs used in our work are crystal structures and each of them has sixteen guanines. Thus, we analyzed six types of dihedrals of sixteen guanine deoxyribonucleotides including experimental structure shown in Fig. 10 and Fig. 11. For potassium-containing GQ, basically all force fields have a tendency to maintain the experimental dihedrals. But some of them deviate largely such as parmbsc0 force field to describing delta and gamma dihedrals. It also signified that the refined parmbsc1 and OL15 are superior to parmbsc0 indeed. In Drude and OL15 force fields, most dihedrals are close to experiment without great deviation. For sodium-containing GQ, AMOEBA and OL15 emerge to deviate a lot towards some dihedrals. However, dihedrals under Drude force field have small deviation to experiment.</p><!><p>Given that QM method is capable of calculating molecular properties authoritatively, it is necessary to contrast the outputs calculated by QM and force field respectively. Thus, we calculated the conformational energy of 20 modified potassium-containing GQs. As shown in Fig. 12, the results of all force field calculations are not strongly correlated with that of QM calculations. AMOEBA force field has the highest correlation, while Drude has the lowest correlation. From above analysis, Drude force field behave much better than other force fields in many aspects. If Drude force field is robust enough to characterize GQ, the conformational energy of GQ obtained by Drude should be congruous with QM calculation as much as possible. It is tricky to explain such a result. However, the outstanding performance of Drude is undeniable as analyzed above.</p><p>Fig. 12 The conformational energy of 20 modified GQs.</p><!><p>We tested and evaluated the performance of five force fields that are parmbsc0, parmbsc1, OL15, Drude2017 and AMOEBA through MD simulation. Drude2017 surpassed other four force fields from multiple aspects except for some imperfection such as conformational energy calculation.</p><p>Currently, none of force fields is absolutely accurate, but Drude2017 has made a great stride to describe GQs. It is a reliable force field to simulate GQ more exactly. It includes induced polarization effect, which is pivotal to describe GQ correctly. For Amber force fields, the lack of polarization effect leads to the formation of bifurcated hydrogen bonds and inability to reserve channel ions in loops. AMOEBA is likewise a polarizable force field, but it cannot sustain GQ structures well despite it characterizes DNA excellently. 43 Thus, we recommend that Drude polarizable force field should be considered first for theoretical calculation of GQs.</p><p>As GQ has more and more extensive application prospects, using solid force fields to conduct MD simulation of GQ is an imperative challenge. In light of outstanding performance of Drude polarizable force field, there is likely an evolution from tuning parameters of pair-additive force fields to developing novel polarizable force field over next few years, in spite of the difficulties such as great computational consuming. More precise and widely accepted polarizable force fields may be present in the future.</p>