C-A test of DNA force fields

The DNA duplex may be locally strongly bent in complexes with proteins, for example, with polymerases or in a nucleosome. At such bends, the DNA helix is locally in the non-canonical forms A (with a narrow major groove and a large amount of north sugars) or C (with a narrow minor groove and a large share of BII phosphates). To model the formation of such complexes by molecular dynamics methods, the force field is required to reproduce these conformational transitions for a naked DNA. We analyzed the available experimental data on the B-C and B-A transitions under the conditions easily implemented in modeling: in an aqueous NaCl solution. We selected six DNA duplexes which conformations at different salt concentrations are known reliably enough. At low salt concentrations, poly(GC) and poly(A) are in the B-form, classical and slightly shifted to the A-form, respectively. The duplexes ATAT and GGTATACC have a strong and salt concentration dependent bias toward the A-form. The polymers poly(AC) and poly(G) take the C- and A-forms, respectively, at high salt concentrations. The reproduction of the behavior of these oligomers can serve as a test for the balance of interactions between the base stacking and the conformational flexibility of the sugar-phosphate backbone in a DNA force field. We tested the AMBER bsc1 and CHARMM36 force fields and their hybrids, and we failed to reproduce the experiment. In all the force fields, the salt concentration dependence is very weak. The known B-philicity of the AMBER force field proved to result from the B-philicity of its excessively strong base stacking. In the CHARMM force field, the B-form is a result of a fragile balance between the A-philic base stacking (especially for G:C pairs) and the C-philic backbone. Finally, we analyzed some recent simulations of the LacI-, SOX-4-, and Sac7d-DNA complex formation in the framework of the AMBER force field.


I. INTRODUCTION
Understanding the mechanisms of the protein binding to DNA is a critical task not only in theoretical biophysics, but also in medicine. To date, there appeared a lot of reports on modeling the processes of DNA-protein binding or dissociation, protein diffusion and conformational dynamics in complexes (see reviews [1,2]). Usually, the researchers choose the AMBER force field. This force field allows simulating the complexes even in the case * zubova@chph.ras.ru of highly deformed DNA (with kinks as a result of protein intercalation), for example, in complexes with Sox protein [3][4][5], Hbb [6,7], Sac7d [8][9][10],LacI [11] and the others [12,13]. The CHARMM force field is being used much less frequently, with the important exception of nucleosomes [14][15][16].
However, the reliability of the results of such modeling is limited by the accuracy of the phenomenological force fields used for protein, DNA and their interactions with each other and with the solvent. Both the AMBER and CHARMM DNA force fields are known to have several drawbacks. Under the conditions when the B-form is to be stable, on time intervals greater than a microsecond, a long DNA molecule loses its duplex structure in the framework of the CHARMM36 force field, and passes into an intermediate form between the B-and A-forms in the framework of the AMBER bsc0 force field [17]. In the 2010s, several torsion angles were edited (primarily ǫ/ζ) to improve the representation of the BII conformation of phosphates [18,19] (AMBER bsc1 and CHARMM36). Recently, in 2021, a correction of the α/γ angles has been proposed to more adequately reproduce the Z-DNA in the AMBER force field (OL21) [20].
More serious shortcomings for practical purposes are too weak base pairing and too strong base stacking in both the AMBER and CHARMM force fields, with the pairing weaker in the CHARMM force field, and the stacking stronger in the AMBER force field (see, for example, the discussion in Ref. [21]). Of the nonbonded interactions, DNA-ions interactions (in particular, with Mg 2+ ) also do not correspond to the experimental ones. The appropriate corrections for the AMBER force field has resulted in a new force field DES-Amber [22] with adjusted values for a number of van der Waals parameters and for some atomic partial charges. In the framework of the DES-Amber force field, one mainly uses a 'native' water model, TIP4P-D, since the TIP3P water is known to be a poor solvent for proteins [23].
In this regard, let us remind that it was the TIP3P water model that was used to derive the values of partial charges on DNA atoms in the CHARMM force field, and these values have never been revised. In the AMBER force field, a new set of partial atomic charges for DNA and RNA has been devised [24]. In order to amend excessively strong electrostatic interactions in DNA-protein complexes, a new Tumuc1 force field has been developed [25], with new partial atomic charges and parameters of the valence bonds, angles and torsion angles. Unexpectedly, the Tumuc1 force field simply adapts all the van der Waals parameters from the AMBER force field.
One can do the opposite and weaken too strong DNAprotein interactions by amending their van der Waals, rather than electrostatic, interactions. The CUFIX (or NBFIX) correction [26] adjusts the σ parameters of the Lennard-Jones potentials between the atoms of the DNA phosphate group and of some groups of atoms on amino acids. The analogous corrections have been offered also for interactions with ions. It turned out that for each pair of atoms, the parameter σ has to be chosen independently of other pairs (combination rules do not apply). The CUFIX correction provides the correct values for the diffusion coefficients of proteins along DNA in both the cases of fast and slow diffusion [7].
However, the most difficult task, currently unsolved, is the adequate modeling of the A-form of DNA. DNA takes the A-form locally in complexes with endonucleases, transcription factors, and polymerases [27,28]. This is a form with a narrow major groove and sugars adopting the north conformation (C3'endo). The geometry of the A-DNA is opposite to that of the C-DNA, in which the minor groove is narrow (see Figure 1). The signs of the Roll and Slide parameters are also opposite in the Aand C-forms. In the C-DNA, a significant fraction (over 40%) of phosphates adopts the BII conformation [29]. DNA in the nucleosome has about 20% of phosphates in the BII conformation (see, for example, the structure 1kx5.pdb [30]).
The A-and C-forms of DNA are competitors under conditions with a low dielectric constant of the medium, or a low water activity, or a small available volume. These are the conditions in hydrophobic protein cavities. In vitro, the competition can be demonstrated in experiments when one increases the alcohol proportion in a dilute DNA aqueous solution [35]. If there is no additional salt in the solution, then, at 78% of ethanol, the natural DNA (calf thymus) takes the A-form. But, with only 0.1M of NaCl salt, the addition of methanol causes the transition to the C-form. In complexes with proteins, DNA can locally take the non-canonical A-and C-forms or be strongly bent [36]. In these cases, it is important that, for conditions corresponding to the experiment, the force field adequately reproduces the transition to non-canonical conformations of the torsion angles on the DNA sugar-phosphate backbone, including the transition to non-canonical deoxyribose conformations (usually north, C3'endo).
Unfortunately, the transition to the A-form is not adequately modeled by both the AMBER and CHARMM force fields. They are considered to be B-and A-philic, correspondingly [17,[37][38][39][40]. Since the balance of interactions in the molecule depends on many factors, changing any of them can easily amend the situation in one particular case. It is more important to find out exactly which parameters of the force fields do not correspond to reality.
The answer to this question requires a direct comparison of the results of simulations and experiments on the B to C and B to A transitions. In an adequate force field, the DNA molecule is to take both the C-and A-forms (depending on the base sequence) under the experimental conditions. For the test, one has to choose the experiments with the simplest possible conditions, without the interaction of the DNA with complex molecules (alcohol; spermine; drugs; polymers, including DNA: in fibers or in crystals). The simplest system for molecular modeling is a DNA aqueous solution with the addition of salt. We have analyzed (Section II) the available experimental data on such transitions for NaCl salt (obtained by Raman and IR spectroscopy [41,42] and circular dichroism [43]).
Based on the analysis, we have formulated the "C-A" test for the balance of interactions between the base stacking and the backbone flexibility in DNA force fields. We have tested (Section III) the AMBER bsc1 and CHARMM36 force fields and their "hybrids": the fields with potentials for the sugar-phosphate backbone from one force field, and for the base stacking interactionsfrom another field. For both the AMBER and CHARMM force fields, we have localized the weak points not allowing one to reproduce the experimentally observed DNA behavior. Finally, in the light of the results obtained, we have reviewed (Section IV) the possible influence of the choice of force field on the results of modeling DNA complexes with three different proteins (LacI [11], Sox-4 [5], and Sac7d [9]), where the DNA is bent towards the major groove and locally takes the A-form. visualization: VMD program [33,34]). For the C-DNA, we show the variant with 100% of BII phosphate conformations. We marked with green and red arrows the local narrowing of the DNA minor and major grooves on the C-like and A-like DNA bends on the nucleosome and in the complex with polymerase, respectively. We also show several more protein complexes with the A-like DNA bends (LacI, SOX-4, and Sac7d).
CD spectrum very close to that of the calf thymus DNA. Poly(AC) also passes into the C-form upon addition of 4.5M NaCl [44]. According to spectral measurements, the polymer poly[d(GC)] 2 (poly(GC)) passes into the Zform [45], but, at a low salt concentration (up to 0.2M), poly(GC) is in the classical B-form. Under the same conditions, the polymer poly(dA)·poly(dT) (poly(A)) also takes a slightly deformed B-form with some share of sugars in the north conformation (C3'endo) [46]. The poly(A) double helix has not been observed in any other form under any conditions. A bias toward the A-form is observed in the (double stranded) oligomers d(ATAT) 2 (ATAT) and d(GGTATACC) 2 (GGTATACC) [47,48] (when they are stable: at low temperatures). Only the polymer poly(dG)·poly(dC) (poly(G)) demonstrates the complete transition to the A-form upon addition of the 4M NaCl salt [45]. Lowering the temperature facilitates this transition. Moreover, even in the almost complete absence of the additional salt (0.03M), 30% of this polymer is in the A-form [45].
Although the picture of the conformational transitions described above seems clear, quantitative estimates can currently be obtained only from IR or Raman spectra. In these spectra, there are definite bands -markers of the A-form [41,42]. Unfortunately, using the markers known for the C-form, one cannot make reliable quantitative estimates. As we will see below, the A markers cannot always provide good accuracy as well. The accuracy can be increased by using 2D IR spectra [49][50][51]. However, such data are not available for the required DNA sequences in the required frequency range.
Certain vibrational modes of the DNA nucleotide change their frequency when the sugar conformation changes. Among these modes are, for example, some deformation modes of the pyrimidine ring on bases. Other modes change their frequency when the geometry of the sugar-phosphate backbone changes. When some of DNA nucleotides are in the A-form, and the rest is in the Bform, one may estimate the fraction of the A-form as the ratio of the integral intensity of the A-form marker band to the sum of the intensities of the A-and B-form bands.
The transition from the B-DNA to the A-DNA is accompanied by the shift in the frequency of the mode OPOsym (symmetrical stretching vibrations in the phosphodiester group (C-O3'-P-O5'-C)) from 790 to 807 cm −1 . For DNA molecules consisting only of G:C pairs (for example, poly(G) and poly(GC)), this band can be used -but only at high salt concentrations -to fairly accurately estimate the degree of transition of the backbone to the A-form (see Table. I and Figure 2). Unfortunately, in this region of the spectrum, in addition to this A marker, there is also a band at a frequency of 780 cm −1 , which corresponds to a certain deformation mode of the cytosine ring. It is possible to separate this band from the OPOsym band only if the salt concentration is high (see Figure 2). In other cases, this marker does not allow quantifying the fraction of the A-form without additional assumptions.
The situation is even worse for oligomers containing A:T pairs. There is a band of thymine similar to that of cytosine at 780 cm −1 , but the frequency of this band shifts after the B-A transition. It is possible that the intensity of this band also changes, as does the intensity of the OPOsym band. Therefore, the optical spectra (both IR and Raman) do not allow quantifying the fraction of the backbone that has passed into the A-form; only qualitative evaluation is possible. Therefore, for the majority of oligomers and polymers in the Table I, we denote the estimates for the OPOsym mode and for the mode of thymine at 777 cm −1 with a question mark. In Figure 2, we give only one of many possible curve fitting of the experimental spectrum by a sum of Gaussians, it cannot be considered reliable.
Another problem for with A:T pairs is the coincidence of the frequencies (642 cm −1 ) of the pyrimidine ring breathing modes for the adenine and the thymine. It is impossible to separate these bands (and estimate their widths and intensities in different DNA forms) from Raman or IR spectra. Therefore, in the table I we present the total estimate for this marker.
The polymer poly(GC) is believed to be in the B-form at low salt concentrations (0-0.2M). In the spectrum of poly(A), there is a noticeable shoulder (at 807 cm −1 ) indicating the presence of north sugars. The spectrum of poly(A) also has a large peak at 820 cm −1 (not marked out in Figure 2). This band with different heights is present on all the Raman spectra from the DNA samples. In the work [46], the authors have even introduced a special designation, a2, for a DNA form with certain (non-canonical) torsion angles. The band at 820 cm −1 was a marker of this new form. Later, the form a2 disappears from the works on IR and Raman spectroscopy of DNA. There is no generally accepted assignment of this band.
The oligomer ATAT exists as a duplex only at low temperatures. From the broadening of the band near 790 cm −1 , one can definitely say that the addition of salt shifts the conformation towards the A-form, but one cannot give reliable quantitative estimates. The same can be said about the oligomer GGTATACC.
The most A-philic polymer is poly(G). Even at the lowest salt concentration (0.03M) and a rather high temperature of 30°C, the band OPOsym has a shoulder corresponding to the A-form. At room temperature and a salt concentration of 0.2 M, the backbone of this polymer is 62% in the A-form.
Let us compare the fractions of north sugar conformations calculated using different markers. The highfrequency thymine deformation mode at 1239 cm −1 gives approximately the same estimate of ∼55% for all the cases presented in Table I. However, judging by the backbone marker OPOsym, the polymer poly(A) is in the (slightly deformed) B-form, and the oligomer ATAT undoubtedly contains a significant fraction of the A-form at the high salt concentration. In the works cited in the Table I, there is only one case when the estimate by this marker differs from ∼55%. This is a crystal GGTATACC in the A-form with t(1239):∼67%. A similar band for the guanine deformation mode (breathing of the pyrimidine ring) at 1318 cm −1 behaves quite differently. For the polymer poly(G), the estimate g(1318) practically repeats the estimate by the backbone marker OPOsym (in the cases when this estimate can be considered reliable). In addition, the estimate by the marker [a+t](642) is always very low. Judging by these results, the use of the markers t(1239) and [a+t](642) cannot be considered justified without further analysis of their validity.
The practical conclusion from the experimental data reviewed above is as follows. Two DNA sequences, poly(AC) and poly(G), behave oppositely in aqueous NaCl solution at a high salt concentration. The first passes into the C-form (the minor groove narrows), the second into the A-form (the major groove narrows). Strictly speaking, one may limit oneself to these two sequences to test a DNA force field for the balance between the flexibility of phosphates and sugar puckers. Other oligomers can serve as reference points. The polymers poly(GC) and poly(A) at low salt concentrations are the variants of the B-form (the classical form and a form with a shift to the A-DNA). For the A-philic oligomers ATAT and GGTATACC, the fraction of the A-form is to increase with the addition of salt.

III. C-DNA AND A-DNA IN AQUEOUS SOLUTION: MODELING
We have simulated all the six DNA sequences discussed in the previous section: poly(GC) and poly(A) (20 bp duplexes); ATAT and GTATACC; poly(AC) and poly(G) (20 bp duplexes) at different salt concentrations (see Table V). We have tested four force fields: AMBER (bsc1), CHARMM (36) and two hybrid fields, in which the sugarphosphate backbone was modeled in the framework of one field, and the base stacking -in the framework of the other force field (for details, see Section VI). In all the cases, we used the same water and ion models (TIP4P/long [53] water, Joung-Cheatham ions [54]).

FIG. 2:
Raman spectra in the region of OPOsym (symmetrical stretching in the phosphodiester group) for some regular DNA polymers and oligomers (extracted from [45,47,52]). A conjectured Gaussian peak fitting is done by us.  Table II and figure 3 show the simulation results for the first four oligomers at a low salt concentration. The first defect of all the tested force fields is that none of them (except for the hybrid 'AMBER -CHARMM's bases' in the case of the oligomer ATAT) showed any dependence on salt concentration for the oligomers ATAT and GGTATACC. This contradicts the experiment (higher salt concentration -larger fraction of the A-form). All the system parameters at the high salt proved out to be the same as at the low salt, so they are not listed in Table II and in Figure 3.
In the AMBER force field, at any salt concentration, II: Conformations of the DNA oligomers for which no dependence on the salt concentration was found in the simulations (except for ATAT for the force field 'AMBER -CHARMM's bases'). Initially, all the oligomers were in the B-form. Unless otherwise indicated, the transition to the final conformation occurred within 5 ns after the system preparation (see Section VI). Temperature: 300K (unless otherwise stated), additional salt: ∼0.15M NaCl, water: TIP4P/long [53]; ions: Joung-Cheatham [54]. For each DNA sequence, the average experimental value of the fraction of BII phosphate conformations in solutions is given (by 31 P NMR chemical shift [55]). C3: north sugar pucker (C3'endo) of thymine (t), cytosine (c), guanine (g) or adenine (a). The units of the parameters Roll and Slide are degrees and angstroms, respectively.
terminal H bonds terminal H bonds 5% BII broken after sev. ns off and on these four oligomers retain the B-form, the oligomer ATAT exhibiting an excess of BII phosphate conformations (at low temperature). In the CHARMM force field, the oligomers GGTATACC and ATAT are in the C-form even at low salt. The bias towards the C-form is the defect of the sugar-phosphate backbone of the CHARMM force field, because this bias is absent in the hybrid force field 'AMBER -CHARMM's bases', but is present in the force field 'CHARMM -AMBER's bases'.
On the other hand, in the CHARMM force field the base stacking is A-philic. Indeed, in the hybrid force field 'AMBER -CHARMM's bases', the amount of north sugars is much higher for all four oligomers and is in the best agreement with the experiment. Moreover, this force field reproduces an increase in the fraction of north sugars with an increase in the salt concentration for the oligomer ATAT. At a salt concentration of 5 M, the replicas (we had four of them) pass from one form to another, being in the B-form more often. So one can say that the result is a mixture of the B-and A-DNA. On average, the fraction of north sugars is 24(±16)% (at 5(±9)% BII phosphate conformations; Roll=5°, Slide=-0.7Å).
For the polymers poly(AC) and poly(G), a salt dependence does take place (see Tables III and IV, and Figure 4). The AMBER force field reproduces the transition of the polymer poly(AC) from the B-to the Cform with the addition of salt, but does not reproduce the bias of the polymer poly(G) to the A-form, nor the complete transition to the A-form at high salt concentration. In the AMBER force field, this polymer retains the B-form at any salt. Moreover, a decrease in temperature increases the bias of poly(G) in the opposite direction: towards the C-form.
In the CHARMM force field, the polymer poly(AC) is never in the B-form, even with zero additional salt. There are successive transitions from B to C and vice versa. At the high salt concentration, a complete transition to the C-form occurs. The polymer poly(G) is completely in the A-form already at a salt of 0.15M (which does not correspond to the experiment). As the salt concentration increases, the fraction of the A-form decreases (also in contrast to the experiment), while the fraction of BII III: Dependence of the conformation of the polymer poly(AC) on NaCl concentration in different force fields (in the experiment, one observes the transition from the B-to the C-form with increasing salt). Simulation parameters and notation are the same as in Table II.  IV: Dependence of the conformation of the polymer poly(G) on NaCl concentration in different force fields (in the experiment, one observes the full transition to the A-form with increasing salt). Simulation parameters and notation are the same as in Table II.  conformations increases.
In the force field 'AMBER -CHARMM's bases', the polymer poly(G) is in the A-form, regardless of the amount of salt, which confirms the bias of the base stacking (rather than the backbone) in the CHARMM force field to the A-form. In the CHARMM force field at the high salt, the narrowing of the minor groove and the partial transition to the B-and then to the C-form (the fraction of phosphates in the BII conformation increases from 3 to 16%) is provided by the sugar-phosphate back-  Table II. The experimental data on the fraction of the A-form are taken from Table I. bone. Indeed, in the force field 'CHARMM -AMBER's bases' (in which the AMBER bases do not tend to A-like base stacking), even at the low salt, the polymer poly(G) is in the C-form. The bias of the DNA backbone in the CHARMM force field toward the C-form may be provided by the interaction of charges on the backbone with ions in solution. At the high salt, a large amount of sodium ions appear in the minor groove of poly(G) near the atom O4' of the sugar ring (direct contact) and near the atom N2 of guanine. Perhaps the reason is that the atom O4' in the CHARMM force field has a charge of -0.5 versus -0.3691 in the AMBER force field. These additional sodium ions in the minor groove, attracting phosphate groups, seem to provide narrowing of the minor groove and the shift in the balance toward the C-form.
In summary, the base stacking of the AMBER force field is excessively B-philic with a slight bias towards the C-form. Contrary to experiment, there is no A-DNA for the polymer poly(G) in this force field. However, it reproduces the transition to the C-form for the polymer poly(AC). In the CHARMM field, the base stacking pro-vides the bias to the A-form for poly(G), while the backbone is excessively C-philic. It can be said that the B-DNA in the CHARMM force field is the result of a fluid balance between the shift of the base stacking towards the A-form and the backbone -towards the C-form. In contrast, the B-DNA in the AMBER force field is a stable inflexible construct due to the B-philic excessively strong base stacking. But the researcher may erroneously conclude that the binding mechanism is a conformational selection (a protein binds only to a certain non-canonical form of DNA), while in reality an induced fit takes place (a protein induces the needed change in DNA conformation during or after binding, see review [1]). Let us compare the results of modeling of the DNA/LacI [11], DNA/SOX-4 [5] and DNA/Sac7d [9] complexes. All the authors use the AMBER force field (variants bsc1 in Ref. [11] and bsc0 in Ref. [5,9]) with the TIP3P water. The bsc1 [18] modification increases the fraction of BII phosphate conformations (underestimated in bsc0) and facilitates the transition to the C-form.
In the experiment, in the process of the complex formation, all three the proteins LacI, SOX-4, and Sac7d bind to the minor groove and, expanding it, form an Alike bend towards the major groove (see Figure 1). DNA locally takes the A-form. In the case of the protein LacI, the binding (and presumably the target recognition) also occurs in the major groove.
In the modeling, the protein Sac7d (albeit at an elevated temperature of 338K) quickly (within tenths of a nanosecond) penetrates the DNA minor groove, moves along it, and induces a bend on the target, showing classical induced fit within a time interval less than a microsecond. The protein SOX-4 also demonstrates fairly fast (within 5 nanoseconds at 300K) non-specific binding. However, for this protein at a lower temperature -and, perhaps, due to too strong DNA-protein interactions in the AMBER force field -the transition from the nonspecific to specific binding seems impossible to observe by ordinary molecular dynamics. The authors do not do this and study the specific binding mechanism with the use of the metadynamics method.
It should be noted that, in most works on DNA-protein  Table I. complex formation, the authors use enhanced sampling techniques: steered molecular dynamics, metadynamics, umbrella sampling, or weighted histogram analysis method. These approaches make it possible to build free energy profiles and compare the obtained energy of complex formation with experiment. But one can not follow the process of complex formation in real time. In case of insufficient number or poor choice of reaction coordinates, one may face a difficulty in finding the real pass to the global minimum. The discrepancies may also result from an inadequate balance of interactions in other parts of the system (for example, incorrect DNA-protein binding energy, inappropriate water model, and so on, see Section I).
In the case of the protein LacI, after a microsecond ordinary MD calculation, the authors failed to obtain the native state of the DNA-protein complex with an A-like bend. Contrary to the proteins SOX-4 and Sac7d, the transition to this state is impeded by a steric clash between LacI and the straight DNA. Because of the known B-philicity of the AMBER force field, the authors hypothesized that the DNA molecule is to be already bent before the complex formation (conformational selection).
However, there are two other factors (specifics of the complex formation and the features of the force field not related to DNA) that could contribute to the negative result of the simulation.
There are experimental data [56] that two "hinge" alpha-helices, bound with the DNA minor groove in the LacI/DNA complex, are disordered in the free and nonspecifically bound protein. More exactly, these data refer to a truncated protein, the "headpiece", consisting of the hinge helices and the binding domains (see Figure 1). But, since the affinity of the headpiece for DNA is the same as of the whole LacI protein, the headpiece is believed to be an adequate model system for binding of the whole protein. In the simulations by Ref. [11], the hinge helices in the free (whole) protein are almost ordered. In the CHARMM force field (also with the TIP3P water), isolated hinge helices in solution are disordered [57]. But they are ordered as a part of the free LacI headpiece [57], contrary to experiment. Such an excessive stability of the hinge helices may result from using the TIP3P water model which is a poor solvent for proteins [23]. The same may be true also for the modeling of the whole LacI protein in Ref. [11].
The ordered state of the hinge helices may not facilitate, but, on the contrary, hinder the formation of the complex. Indeed, the folded hinge helices cannot penetrate the minor groove at the early stages of binding. So they cannot destroy the spine of hydration [58] in the minor groove, which stabilizes the B-form and prevents binding as well as bending. Without additional research, it is impossible to say what exactly -B-philic DNA of the AMBER force field, inappropriate water model, both together, or something else in addition (as a high barrier for the binding) -caused the negative result of modeling in Ref. [11].
In the cases of SOX-4 and Sac7d, such a complicated rearrangement of the protein is not needed for the complex formation. Therefore, excessively strong DNAprotein interactions may outweigh the DNA rigidity in the AMBER force field. As a result, the protein Sac7d easily induces the required DNA bend. And, for the native DNA/SOX-4 complex, there is a corresponding (global) minimum and a not very high barrier to pass from the (local) minimum for the non-specifically bound complex.
Thus, when analyzing the results of modeling the formation of DNA-protein complexes, it is necessary to take into account both the specifics of the possible mechanism of binding of the particular protein to DNA, and the features of the force fields used. Unfortunately, the choice of the force field is almost never substantiated in the corresponding works. In the case of modeling the complex with the protein LacI, it would be better to choose the CHARMM force field for the DNA. In it, the DNA duplex is more flexible and, for sequences with a large number of G:C pairs, more A-philic. This choice would allow eliminating at least one of the uncertainty factors in the analysis of the obtained simulation results.

V. CONCLUSION
The DNA duplex in the DNA-protein formation process often locally bends (nucleosome, complexes with polymerases) and passes into the C (narrow minor groove, negative Roll, positive Slide) or A (narrow major groove, positive Roll, negative Slide) forms. The B-C transition is accompanied by the switch of a significant fraction (more than ∼40%) of phosphates from the BI to the BII conformation, and the B-A transition is accompanied by the switch of sugars from the south (C2'endo) to the north (C3'endo) conformation. To study the formation of protein complexes with bent DNA by molecular dynamics modeling, the DNA force field is needed to adequately reproduce both the conformational transitions, which requires appropriate testing.
In the present work, we have analyzed the available experimental data on the B-C and B-A transitions under the conditions most conveniently implemented in molecular dynamics modeling: in an aqueous NaCl solution. Unfortunately, from the currently available spectroscopic data (Raman and IR spectroscopy), one cannot extract in most cases reliable quantitative estimates of the fraction of north sugars or of the fraction of the DNA backbone in the A conformation. For the B-C transition, only the circular dichroism data are reliable. And it is also currently impossible to extract the exact fraction of the BII conformations of phosphates from them. However, there are DNA sequences that can be used to test DNA force fields.
We have selected six DNA oligomers that respond differently to an increase in salt concentration. At a low salt concentration, two polymers, poly(GC) and poly(A), take the B-form, classical and slightly deformed with some fraction of north sugars, respectively. The oligomers ATAT (at a very low temperature, when it is stable) and GGTATACC (also at low temperature) have an appreciable bias towards the A-form. The fraction of the A-form increases with increasing salt concentration. At a high salt concentration of 4-5M, the polymers poly(AC) and poly(G) take the C-and A-forms, respectively.
Modeling these six DNA duplexes at the two (low and high) salt concentrations is a good test for the balance between the base stacking and the backbone (phosphate groups and deoxyriboses) conformational mobility in DNA force fields. We have tested the two most commonly used force fields AMBER bsc1 and CHARMM36. To find out how exactly the balance is upset in each case, we also had to test the hybrid force fields, with backbone interactions taken from one of the force fields, and the base stacking -from another.
We showed that the "B-philicity" of the AMBER force field is a consequence of the "B-philicity" of its base stacking. The stacking is excessively strong, and it does not allow the polymer poly(G) to pass to the A-form. The transition to the C-form for the polymer poly(AC) is possible in the framework of the AMBER force field. In the CHARMM force field, the balance between the interactions in the DNA duplex is fundamentally different. Namely, the B-form is the result of an fragile balance between the A-philic base stacking (especially for G:C base pairs) and the C-philic backbone. This equilibrium is strongly shifted towards the C-form with the predominance of A:T pairs in the DNA sequence, and towards the A-form with the predominance of G:C pairs. Therefore, in the CHARMM force field, the polymer poly(G) is always in the A-form, and the fraction of sugars in the north conformation for poly(GC) and poly(A) approximately coincide (in the experiment, poly(A) is more Aphilic). The question on the balance between the base stacking and the backbone conformational flexibility in the real DNA duplex remains open.
In both the force fields, the interactions between the DNA and the ions make much smaller contribution to the overall balance of interactions in the system as compared to the experiment. In simulation of the oligomers ATAT and GGTATACC the experimental dependence on salt concentration is not observed. The influence of ions in both the force fields is especially weak under the conditions of the B-A transition. The fraction of north sugars in poly(G) does not depend on the salt concentration in the framework of neither force field (except for the wrong dependence in the CHARMM force field); while the B-C transition for poly(AC) is reproduced in the AMBER force field, and almost reproduced in the CHARMM force field. Due to the "dynamic" nature of the B-form in the CHARMM force field, at low salt concentrations the polymer poly(AC) incessantly passes between the B-and C-forms, the C-form having too large percentage of BII phosphate conformations.
Thus, we have revealed a discrepancy in the DNA conformational flexibility (transitions to the C-and Aforms) between modeling in the framework of the AM-BER and CHARMM force fields and the experimental data. A researcher choosing a force field might want to take into account the revealed features of the DNA duplex in these two commonly used force fields, along with the specifics of the possible mechanism of the complex formation, especially if the formation requires A-like bending of the DNA duplex. For example, in the case of the LacI/DNA complex, the CHARMM force field would be a better choice than the AMBER force field. In the CHARMM force field, the DNA duplex is more flexible and, with the predominance of G:C pairs, more A-philic. However, the reliability of the results can only be ensured by a force field that can withstand the C-A test.

VI. METHODS: DETAILS OF MD CALCULATIONS
We simulated DNA oligomers in water with ions using the LAMMPS [59, 60] package in the framework of the force fields AMBER bsc1 [18] and CHARMM36 [19] and their hybrids, four force fields in total. In the hybrid force fields, the bonded interactions belonged to the first force field (for example, AMBER in the force field 'AMBER -CHARMM's bases'). The non-bonded interactions in the sugar-phosphate backbone belonged to the same force field, while the van der Waals parameters and partial electric charges on the bases were replaced by their counterparts from the other force field (CHARMM in the force field 'AMBER -CHARMM's bases'). Thus, the base stacking (van der Waals and Coulomb interactions between bases) was modeled in the framework of the other force field. We have compensated the difference in total charges on the bases (they are neutral in the CHARMM and negative in the AMBER) by a change in the charge on the hydrogen atom H1'. In such a splicing of the force fields, the effective potential of base rotation around the glycosidic bond χ turns out to be hybrid. Indeed, the potential depends on hybrid van der Waals and Coulomb interactions between the atoms on the sugar ring (with the parameters of the first force field) and on the bases (with the parameters of the second force field). We have plotted potential energy curves for rotation around the glycosidic bond χ for all four bases bound to the deoxyribose locked in the south conformation in a vacuum. It proved out that these curves lie between the corresponding curves for the CHARMM and AMBER force fields in the area of the anti− conformation. For calculation of the electrostatic interactions, we used the PPPM (particle-particle particle-mesh Ewald) [61] method with an accuracy of six significant digits (the dimensionless LAMMPS parameter 10 −5 ). We chose to cut off the short-range part of the potential at the distance optimal for the calculation speed (8Å). In the AM-BER force field, we calculated the van der Waals forces with an accuracy of 10 −5 kcal/mol/Å (i.e., we considered the force equal to zero at a distance larger than 8-20Å depending on the type of atoms). In the CHARMM force field, we applied this approach only to ions and water (to their interactions with each other and with the DNA). The van der Waals interactions between the DNA atoms smoothly vanished in the range from 12 to 13Å, regardless of the type of atoms.
We used the Langevin thermostat and the Nose-Hoover barostat (pressure of 1 atmosphere) with a damping parameter of one picosecond. We restrained X-H valence bonds and H-X-H valence angles (X is any atom) with 10 −5 kcal/mol tolerance by the SHAKE algorithm [62]. The integration step was equal to 2 fs.
We set up the initial state of the system as follows. We generated the DNA molecule in the B-form by the code of the D.Case group [63]. We placed sodium and chloride ions at random positions around the DNA no closer than 3Å to any other ion or DNA atom. Randomly oriented water molecules were placed in positions with oxygen atom not closer than 2.5Å to any other atom in the system. Then we corrected the atom positions by energy minimization using the conjugate-gradient method (until the force vector norm became less than 10 3 kcal/mol/Å). After that, we launched one MD run at a temperature T=300K for 10 picoseconds in the NVT ensemble and then another run at the same temperature for 10 picoseconds in the NPT ensemble. One could start the MD relaxation of the system from this state, but the cloud of ions usually comes to equilibrium much later than other parts of such systems. In order to speed up this process, we performed an additional MD run with the fixed DNA molecule in the NPT ensemble for 2 nanoseconds. For the first nanosecond, the thermostat temperature for water molecules was 300K, and for ions -500K. During the next nanosecond, we gradually cooled the ions (and the water if needed) down to the simulation temperature (300K, 283.15K, or 271.15K). The productive runs were carried out in the NPT ensemble.
In long-term all atom simulations, one observes too often and too long openings of the terminal base pairs [64,65]. To avoid the fraying, we linked the bases in the terminal base pairs by a parabolic potential. More exactly, we replaced the energy of one of the hydrogen bonds of the pair (the central one: N1-H3 on A:T and N3-H1 on G:C) with a parabola. In the region where the energy reaches its minimum (2.2 -2.3Å), the parabola approximates the sum of the van der Waals and Coulomb interactions between the nitrogen and the hydrogen atoms. For long oligomers, such restraints on the terminal base pairs should have little effect on the result of simulations.
We simulated the tetramer ATAT without the restraints on the terminal base pairs. It is known that at room temperature this oligomer is unstable and decomposes into two strands. At the experimental temperature of -2°C, when the oligomer is intact, in the force fields AMBER and 'AMBER-CHARMM's bases', the hydrogen bonds did not break during the simulation (55-219 ns in total from 4 replicas). In the CHARMM force field, the terminal hydrogen bonds broke after a few nanoseconds. In the case of the force field 'CHARMM-AMBER's bases' (with stronger AMBER's base stacking), the breaks occurred less frequently and healed in most cases. However, in both the cases, we calculated the average values of the parameters Roll and Slide only on the central step of the oligomer ATAT.
In all the DNA force fields we used the same TIP4P/long water [53] and the same Joung-Cheatham ions [54]. Due to the specifics of the choice of partial charges in the CHARMM force field, the water model could have a significant impact on the calculation results. Therefore, we carried out several additional simulations with the inherent models for ions and water (file toppar water ions.str in CHARMM36). We simulated the oligomer GGTATACC with restrained hydrogen bonds on the terminal base pairs at two salt concentrations 0.15 and 4-5M. The replacement of the water and ions led to a decrease in the fraction of BII phosphate conformations, but only by 20% from ∼60 to ∼40%. The fraction of the north sugars increased from ∼3 to ∼15% at the high salt concentration and up to 19% at the low salt concentration. Both the changes were in the correct direction, but the observed values are still far from the experiment (17% BII and ∼(20-50)% C3'endo at 5-6M NaCl).
The computational cell was cubic, with periodic boundary conditions. The edge of the cell was about 15Å greater than the axis of the extended oligomer in the Bform. We carried out the calculations for the oligomers with counterions, and in weak and strong salt solutions, with concentrations of approximately 0.1-0.15M and 4-5M. Table V lists the total simulation time and the characteristics of the modeled systems depending on the salt concentration. In some cases, we collected statistics from several (usually 4) replicas. We stopped the calculation after the transition to another conformation (if any) was completed, and the lifetime in the new conformation was long enough to calculate the equilibrium average values and their dispersions for the parameters of the DNA helix.