Nohad Gresh
email: gresh@lct.jussieu.fr

Sehr Naseem-Khan

Louis Lagardère

Jean-Philip Piquemal

Judit E Sponer
email: sponer@ncbr.muni.cz

Jiri Sponer

Judit E Sponer

Jiri Sponer

Channeling Through Two Stacked Guanine Quartets of One and Two Alkali Cations in the

come   L'archive ouverte pluridisciplinaire

second-order (polarization, charge-transfer and dispersion) contributions within ∆E is mandatory. With two cations in the channel, the relative weights of the second-order contributions increase steadily upon increasing the ion size. In the G8 complexes with two K + or two Rb + cations, the sum of polarization and charge-transfer exceeds the first-order along terms for all ion positions.

Introduction.

The present study bears on the channeling of one and two alkali cations between two stacked guanine quartets (GQ). It is motivated by the growing evidence for the involvement of G quartets in a wealth of biological events [for a recent review, 1]. Thus GQs promote the formation of quadruple helices in telomeric DNA [START_REF] Williamson | G-Quartet Structures in Telomeric DNA[END_REF] and in G-rich sequences in the genome [START_REF] Rodriguez | Small Molecule-Induced DNA Damage Identifies Alternative DNA Structures in Human Genes[END_REF][START_REF] Biffi | Quantitative Visualization of DNA G-Quadruplex Structures in Human Cells[END_REF] Evidence is mounting for the involvement in GQs in human genome regions that are essential for replication, such as c-myc [START_REF] Siddiqui-Jain | Direct Evidence for a G-Quadruplex in a Promoter Region and its Targeting with a Small Molecule to Repress c-Myc Transcription[END_REF] , bcl-2 [START_REF] Dexheimer | Deconvoluting the Structural and Drug-Recognition Complexity of the G-Quadruplex-Forming Region Upstream of the bcl-2 P1 Promoter[END_REF] , c-kit [START_REF] Wei | Crystal Structure of a c-Kit Promoter Quadruplex Reveals the Structural Role of Metal Ions and Water Molecules in Maintaining Loop Conformation[END_REF] , VEGF [START_REF] Agrawal | Solution Structure of the Major G-Quadruplex Formed in the Human VEGF Promoter in K(+): Insights into Loop Interactions of the Parallel G-Quadruplexes[END_REF] , and Refs. therein] but also in viruses [START_REF] Métifiot | G-Quadruplexes in Viruses: Function and Potential Therapeutic Applications[END_REF][START_REF] Amrane | Topology of a DNA G-Quadruplex Structure Formed in the HIV-1 Promoter: A Potential Target for Anti-HIV Drug Development[END_REF] . GQs thus constitute an emerging target for the design of novel chemotherapeutic molecules which can stabilize the quadruple helix and inhibit telomerase or interfere with transcription 11-15 . GQs are endowed with electron transport [START_REF] Livshits | Long-Range Charge Transport in Single G-Quadruplex DNA Molecules[END_REF][START_REF] Ma'ani Hessari | Programmed Self-Assembly of a Quadruplex DNA Nanowire[END_REF] , photoelectronic properties of potential interest in nanochemistry [START_REF] Jissy | Nucleic Acid G-Quartets: Insights into Diverse Patterns and Optical Properties[END_REF][START_REF] Hua | Cation Effect on the Electronic Excited States of Guanine Nanostructures Studied by Time-Resolved Fluorescence Spectroscopy[END_REF] while GQ-based nanowires are endowed with ion-conducting properties [START_REF] Liu | Direct Measurement of Electrical Transport Through G-Quadruplex DNA with Mechanically Controllable Break Junction Electrodes[END_REF][START_REF] Trajkovski | Assessing Roles of Cations in G-Quadruplex-Based Nanowires by NMR[END_REF] .

GQ quadruple helices have been investigated by X-ray crystallography [START_REF] Parkinson | Crystal Structure of Parallel Quadruplexes from Human Telomeric DNA[END_REF][START_REF] Campbell | Structural Basis of DNA Quadruplex Recognition by an Acridine Drug[END_REF] and Refs. therein, high-resolution NMR [START_REF] Bouaziz | A K Cation-Induced Conformational Switch within a Loop Spanning Segment of a DNA Quadruplex Containing G-G-G-C Repeats[END_REF][START_REF] Lim | Structure of the Human Telomere in K(+) Solution: a Stable Basket-Type G-Quadruplex with only two G-Tetrad Layers[END_REF][START_REF] Agrawal | The Major G-Quadruplex Formed in the Human BCL-2 Proximal Promoter Adopts a Parallel Structure with a 13-nt Loop in K(+) Solution[END_REF] , and molecular dynamics [START_REF] Ren | Tiny Telomere DNA[END_REF][START_REF] Fadrná | Single Stranded Loops of Quadruplex DNA As Key Benchmark for Testing Nucleic Acids Force Fields[END_REF][START_REF] Islam | Extended Molecular Dynamics of a c-kit Promoter Quadruplex[END_REF][START_REF] Stadlbauer | Coarse-Grained Simulations Complemented by Atomistic Molecular Dynamics Provide New Insights into Folding and Unfolding of Human Telomeric G-Quadruplexes[END_REF] ; for a review, see Ref. 31. Gquartets have been studied in several high-level ab initio quantum chemistry (QC) calculations bearing on their electronic and cation-binding properties [START_REF] Bazzi | Designing a New Class of Bases for Nucleic Acid Quadruplexes and Quadruplex-Active Ligands[END_REF][START_REF] Zaccaria | The Role of Alkali Metal Cations in the Stabilization of Guanine Quadruplexes: why K+ is the Best[END_REF] , the ranking of competing conformers [START_REF] Šponer | Relative Stability of Different DNA Guanine Quadruplex Stem Topologies Derived Using Large-Scale Quantum-Chemical Computations[END_REF] and the cooperativity of their self-associations [START_REF] Fonseca-Guerra | Telomere Structure and Stability: Covalency in Hydrogen Bonds, Not Resonance Assistance, Causes Cooperativity in Guanine Quartets[END_REF][START_REF] Paragi | Cooperativity in the Self-Assembly of the Guanine Nucleobase into Quartet and Ribbon Structures on Surfaces[END_REF] . Extensions of QC analyses to novel halogenated derivatives of deazaguanine and their complexes with Na + and K + were recently reported [START_REF] Bazzi | Designing a New Class of Bases for Nucleic Acid Quadruplexes and Quadruplex-Active Ligands[END_REF] .

The alkali cations Na + and K + can act as important contributors to the stability of stacked G quartets, owing to their propensity to bind in-between two successive GQs. GQ-containing quadruple helices can interconvert between several alternative forms with short-lived intermediates, the detection of which can evade direct experimental measurements. Longduration molecular dynamics could be ideally suited for such explorations and have already provided insight to unravel candidate interconversion pathways [START_REF] Stadlbauer | Coarse-Grained Simulations Complemented by Atomistic Molecular Dynamics Provide New Insights into Folding and Unfolding of Human Telomeric G-Quadruplexes[END_REF] . There are, however acknowledged limitations to 'classical' molecular mechanics (MM) force-fields for such applications, which could be magnified in the presence of metal cations. Thus a recent study investigated the channeling of a monovalent cation along the Z axis of two and three stacked G quartets, in the presence of another cation fixed along this axis at 1.1 Å above the topmost quartet [START_REF] Gkionis | Ion Binding to Quadruplex DNA Stems. Comparison of MM and QM Descriptions Reveals Sizable Polarization Effects Not Included in Contemporary Simulations[END_REF] . The MM energy profile differed significantly from the DFT-D3 one upon moving the cation closer to the fixed one. This was ascribed to the lack of an explicit polarization contribution (E pol ): a raise of the inter-cation electrostatic repulsion ought to be, to a significant extent, compensated for by a concomitant increase of E pol , since the polarizing field exerted by both cations on the G quartets increases as the two cationic charges get closer together, yet this contribution is absent from classical MM. This then raises an essential issue: how well, by contrast, would anisotropic polarizable molecular mechanics (APMM) profiles fare with respect to the QC profiles? Owing to the highly polar and polarizable nature of guanine, cation-bound G quartets constitute a very revealing benchmark for the accuracy of APMM procedures. Their stabilization involves the interplay of multiply H-bonded, stacked, and cation-ligand interaction and the onset of non-additivity of the total intermolecular interaction energies (∆E). Evaluation of the APMM accuracy in light of these tests is also part of the construction of polarizable potentials for DNA and RNA [START_REF] Gresh | Theoretical Study of Binding of Hydrated Zn(II) and Mg(II) Cations to 5'-Guanosine Monophosphate. Toward Polarizable Molecular Mechanics for DNA and RNA[END_REF] owing to the prevalence and constant interplay of such effects. The present GQ tests will bear on the SIBFA APMM procedure [START_REF] Gresh | Anisotropic, Polarizable Molecular Mechanics Studies of Inter-and Intramolecular Interactions and Ligand-Macromolecule Complexes. A Bottom-Up Strategy[END_REF] and is a continuation of two preceding works. The first bore on the calibration of the alkali cation series (Li + , Na + , K + , Rb + , Cs + ) and QC validation tests of a series of their polyligated complexes [START_REF] Dudev | Quantum-Chemistry Based Calibration of the Alkali Metal Cation Series (Li+-Cs+) for Large-Scale Polarizable Molecular Mechanics/Dynamics Simulations[END_REF] . The second bore on several stacked and H-bonded complexes of cytosine and guanine [START_REF] Gresh | Analysis of the Intermolecular Interaction Energies by Parallel Quantum Chemistry and Polarizable Molecular Mechanics[END_REF] . The first polarizable molecular mechanics study bearing on the selective entrapment of alkali cations between two G quartets was a SIBFA study published in 1986 [START_REF] Gresh | A Theoretical study of the Selective Entrapment of Alkali and Ammonium Cations Between Guanine Tetramers[END_REF] . Each quartet was then however considered as one entity, while in the present study the individual entities are the bases and the cations.

A CPK representation of the two stacked G quartets and of the channeling cation in its successive positions are represented in Figure 1. The unmoving cation is on the top of the Z axis. Throughout this study we retain for both QC and SIBFA computations the idealized stacked GQ geometry constructed in Ref. [37]. Owing to the validation purpose of the present study, relaxing the geometries by energy-minimization or molecular dynamics would be unrealistic and untimely, as it would ignore the averaging effect of additional stacked G quartets, of other channeling cations and the effects of the environments: this will be very continuation of the present work. Note that the geometries taken from ref. [37] are very well prepared and should be sufficient for the purpose of the benchmark computations.

The organization of this paper is as follows. We first consider one G quartet and compare SIBFA and QC results at both HF and DFT-D3 levels. GQ is emblematic of cyclic, multiply H-bonded complexes with strong anticipated cooperativity. Energy Decomposition Analyses (EDA) performed at the HF level will enable to compare the summed values of each individual energy contribution in the six guanine dimers to its value in the quartet. We thus trace back the share of each contribution of ∆E(QC) in the total cooperativity, and how well it is reflected in the context of the SIBFA. Another test bears on the stacking of two halfquartets. We retain in the first quartet one H-bonded guanine dimer, denoted as G1-G2, and in the second quartet one of the two H-bonded dimers with maximal overlap with G1-G2, denoted as G1'-G2'. This will enable to evaluate non-additivity in mixed H-bonded/stacked complexes: could cooperativity now possibly revert into anticooperativity? Such analyses are limited to four guanines as larger-sized complexes are not amenable to EDA with large basis sets.

We next consider the complex of two stacked G quartets. It involves two extensive surfaces, each guanine of the first quartet being involved in non-remote stacking interactions with two guanines of the second quartet, whence a total of eight such interactions. Two stacked G quartets total twenty-eight simultaneous H-bonded or stacking guanine-guanine interactions, whether close or remote. Could an acceptable match of ∆E(SIBFA) to ∆E(QC) still be obtained? Following these preliminary validations, we consider the channeling of one alkali cation, Li + , Na + , K + , and Rb + along the Z axis. We follow the same trajectory as in the original paper [START_REF] Gkionis | Ion Binding to Quadruplex DNA Stems. Comparison of MM and QM Descriptions Reveals Sizable Polarization Effects Not Included in Contemporary Simulations[END_REF] . It starts from a Z= -1.1 Å position below the first quartet and ends at Z= 2.5 Å above it, namely 0.8 Å below the second quartet. The same trajectory is recomputed, now in the presence of a second cation, the position of which is fixed at 1.1 Å above the second quartet. The shortest M + -M + distance of approach is thus 1.9 Å. The first series of SIBFA vs. QC comparisons are reported at the HF level and with HF-derived multipoles and polarizabilities for SIBFA calculations. The second series of comparisons is at the DFT-D3 level and, for the SIBFA calculations, with DFT-derived multipoles and polarizabilities and the inclusion of an explicit dispersion contribution. We conclude with perspectives of MD applications to GQ-based channels and their alkali cation complexes.

Procedure.

Quantum-Chemistry (QC) calculations. The QC calculations used the cc-pVTZ(-f) basis set 43, 44 since the studied complexes on account of their sizes were not amenable to aug-cc-pVTZ(-f) calculations. A triple-dzeta basis set was used on K + 45 . The Stuttgart effective core potential [START_REF] Bergner | Ab Initio Energy-Adjusted Pseudopotentials for Elements of Groups 13-17[END_REF] was used on Rb + . Energy Decomposition Analyses (EDA) resorted to the RVS analysis due to Stevens and Fink [START_REF] Stevens | Frozen Fragment Reduced Variational Space Analysis of Hydrogen Bonding Interactions. Application to the Water Dimer[END_REF] and coded in the GAMESS-US package [START_REF] Schmidt | General Atomic and Molecular Electronic Structure System[END_REF] . Calculations at the correlated level used the B97-D3 energy functional augmented with the dispersion correction by Grimme et al. [START_REF] Grimme | A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu[END_REF]. Some DFT calculations were also performed with the B3LYP functional [START_REF] Lee | Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density[END_REF][START_REF] Becke | A New Mixing of Hartree-Fock and Local Density-Functional Theories[END_REF] . The DFT calculations were done with the G09 package [START_REF] Frisch | Revision B.01[END_REF] .

SIBFA calculations.

In the SIBFA procedure 39, the intermolecular interaction energy is computed as the sum of five contributions: electrostatic multipolar (E MTP* ), short-range repulsion (E rep ), polarization (E pol ), charge transfer (E ct ), and dispersion (E disp )

∆E TOT = E MTP* + E rep + E pol + E ct + E disp
E MTP* is computed with distributed multipoles (up to quadrupoles) derived from the QC molecular orbitals precomputed for each individual molecule using the Generalized Multipole Analysis (GDMA) method by Stone [START_REF] Stone | Distributed Multipole Analysis: Stability for Large Basis Sets[END_REF] . It is augmented with a penetration term [START_REF] Piquemal | Improved Formulas for the Calculation of the Electrostatic Contribution to the Intermolecular Interaction Energy from Multipolar Expansion of the Electronic Distribution[END_REF][START_REF] Piquemal | Toward a Separate Reproduction of the Contributions to the Hartree-Fock and DFT Intermolecular Interaction Energies by Polarizable Molecular Mechanics with the SIBFA Potential[END_REF] . The anisotropic polarizabilities intervening in the expression of E pol are distributed on the centroids of the localized orbitals (heteroatom lone pairs and bond barycenters) using a procedure due to Garmer and Stevens [START_REF] Garmer | Transferability of Molecular Distributed Polarizabilities from a Simple Localized Orbital Based Method[END_REF] . E rep and E ct , the two short-range repulsions, are computed using representations of the molecular orbitals localized on the chemical bonds and on localized lone-pairs [START_REF] Piquemal | Toward a Separate Reproduction of the Contributions to the Hartree-Fock and DFT Intermolecular Interaction Energies by Polarizable Molecular Mechanics with the SIBFA Potential[END_REF][START_REF] Gresh | Model, Multiply Hydrogen-Bonded Water Oligomers (N = 3-20). How Closely Can a Separable, ab Initio-Grounded Molecular Mechanics Procedure Reproduce the Results of Supermolecule Quantum Chemical Computations?[END_REF] . E disp is computed as an expansion into 1/R**6, 1/R**8, and 1/R**10, and also embodies an explicit exchange-dispersion term [START_REF] Creuzet | Adjustment of the SIBFA Method for Potential Maps to Study Hydrogen Bonding Vibrational Frequencies[END_REF] .

We have also calibrated the short-range cation-cation repulsion, as it was found that at shorter M + -M + distances, the sole electrostatic repulsion was insufficient. The calibration bore on the atom-pair multiplicative factor of E rep , denoted PK IJ , where I and J denote the atomic numbers of the two atoms concerned. It was done on the basis of EDA computations on the M+-M+ diatomic complexes, to match the radial dependency of E rep in a range of distances between 2.0 and 3.0 Å. The multiplicative factors are PK Li-Li =1.0, PK Na-Na =9.0, PK K-K =300.0, and PK Rb-Rb =685. Apart from this point, we use the same cation parameters as those derived in 2015 [START_REF] Dudev | Quantum-Chemistry Based Calibration of the Alkali Metal Cation Series (Li+-Cs+) for Large-Scale Polarizable Molecular Mechanics/Dynamics Simulations[END_REF] .

The multipoles and polarizabilities at the HF and DFT levels are those derived in the context of our above-mentioned study on the cytosine and guanine complexes [START_REF] Gresh | Analysis of the Intermolecular Interaction Energies by Parallel Quantum Chemistry and Polarizable Molecular Mechanics[END_REF] . The general parameters were those previously derived on the basis of EDA on model complexes with the aug-cc-pVTZ(-f) basis by Devereux et al. [START_REF] Devereux | A Supervised Fitting Approach to Force Field Parametrization with Application to the SIBFA Polarizable Force Field[END_REF] with the automated I-NoLLS (Interactive, nonlinear least squares) procedure [START_REF] Law | I-NoLLS: A Program for Interactive Nonlinear Least-Squares Fitting of the Parameters of Physical Models[END_REF] . Following [El Khoury et al., submitted] we use an improved representation of the sp 2 lone pairs on O6 and of the sp lone pairs on N7 and N3, 'smearing' them partially on both sides of the guanine plane. The in-plane component of the N7 and N3 lone pairs retains an electron population of 1.45, while each out-of-plane component, extending the external bisector of the nitrogen perpendicular to the plane, has an electron population of 0.275. The in-plane component of each O6 sp 2 lone-pair retains an electron population of 1.70, and there are two smeared components each with an electron population of 0.30, above and below the C=O bond. Such populations, together with the internal coordinates of the smeared lone-pairs, were fit to improve the radial dependencies of E rep (SIBFA) compared to E exch (RVS) upon in-and out-of-plane variations of the approach of a Zn(II) probe to O6 and N7 [Naseem-Khan et al., to be published]. We denote as E 1 the sum of the two first-order contributions E C and E exch (QC) and E MTP and E rep (SIBFA), and by E 2 the sum of the two second-order contributions E pol and E ct .

Non-Covalent Interaction (NCI) analyses. The NCI procedure [START_REF] Johnson | Revealing Non-Covalent Interactions[END_REF][START_REF] Contreras-García | NCIPLOT: a program for plotting non-covalent interaction regions[END_REF] is based on the study of the reduced density gradient as a function of density. It enables to visualize interaction zones between two molecules and gives access qualitatively to the magnitude of these interactions. HF results. The RVS procedure enables to trace back the origin of cooperativity by comparing the value of each ∆E contribution to the sum of its values in the six guanine dimers considered separately. The same comparisons are reported for SIBFA. ∆E(HF) is overestimated by ∆E(SIBFA) by 4.5 kcal/mol out of 50. This is due to both E MTP and E rep within E 1 , each accountable for 2.1 kcal/mol out of app. 65 and 40 respectively. These overestimations occur in the four H-bonded dimers, amounting to 0.6 kcal/mol out of 16 and 10 for E MTP* and E rep , respectively. These errors could appear of lesser importance, but each is multiplied by four upon passing to the quartet. Shortcomings of E MTP* could be due to instabilities in the GDMA procedure to derive distributed multipoles in the case of highly conjugated anionic ligands with diffuse functions, rendering them very sensitive to the choice of input parameters to derive them from the molecular orbital, namely the effective H radii and the switch function [START_REF] Stone | Distributed Multipole Analysis: Stability for Large Basis Sets[END_REF] . Towards deriving less sensitive multipoles, newly-emerging alternatives to GDMA are being considered, such as the Iterative Stockholder Analysis (ISA) 63, 64 and Gaussian Electrostatic Multipoles (GEM) [START_REF] Elking | Gaussian Multipole Model (GMM)[END_REF][START_REF] Cisneros | Application of Gaussian Electrostatic Model (GEM) Distributed Multipoles in the AMOEBA Force Field[END_REF] . Their use will be reported elsewhere. E rep is in turn impacted by the multipolar distribution as it has a prefactor which depends upon the monopolar charges of the interacting atoms. At present we move forward and leave aside these shortcomings. As shown below, upon passing to correlated levels, correlated multipoles afford a much closer match of SIBFA to QC than at the HF level. Furthermore, E 2 (SIBFA) gives a closer match to its RVS counterpart than E 1 , and it is E 2 that is responsible for the non-additivity of ∆E, δE nadd , a critical element in the stacked G quartet-cation complexes.

What is the amplitude of δE nadd (QC), and could it be satisfactorily accounted for? Table 1 also reports the values of ∆E and its contributions in each of the six guanine-guanine dimer complexes, denoted G1-G2 to G3-G4. There are four equivalent doubly H-bonded dimers, G1-G2, G1-G4, G2-G3, G3-G4, and two non H-bonded dimers, G1-G3 and G2-G4. The summed values correspond to four times the values of one H-bonded dimer plus twice that of a non H-bonded dimer. As in our previous papers, we compare E pol (RVS) to E pol *(SIBFA), as computed prior to iterating on the induced dipoles, and E pol (VR) to E pol (SIBFA), after completion of the iterations. E pol (VR) is obtained as the difference between ∆E(HF) and E 1 (RVS) and E ct (RVS). δE nadd (SIBFA) of -9.2 kcal/mol is close to the δE nadd (QC) of -9.6 kcal/mol. A close agreement is also found regarding polarization and charge-transfer. Thus for E pol (RVS) and E pol *, δE nadd amounts to -5.9 and -5.8 kcal/mol (respectively). It increases in magnitude for both E pol (VR) and E pol (SIBFA), to -8.8 and -8.4 kcal./mol, respectively. Such increases were previously noted as signatures for cooperativity. E ct (RVS) is also clearly cooperative, even though its δE nadd has a lesser magnitude (-0.8 kcal/mol in both QC and SIBFA). Such a result has precedents with the cyclic water pentamers and hexamers [57]. In these complexes, as with GQ, each monomer simultaneously acts as proton acceptor from one neighbor and a proton donor to the other. 2). ∆E tot (SIBFA) amounting to -68.8 kcal/mol has a much closer agreement with ∆E(B97-D3) amounting to -69.8 kcal/mol than with ∆E(B3LYP-D3) amounting to -74.7 kcal/mol. Such closer agreements with the former functional have been previously noted [START_REF] Zhang | Polarizable Multipole-Based Force Field for Aromatic Molecules and Nucleobases[END_REF][START_REF] Lemkul | Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model[END_REF] . δE nadd has similar values with both functionals (-11.5 and -11.8 kcal/mol), however δE nadd (SIBFA) is smaller (-8.8 kcal/mol). We have no explanation for this reduction in magnitude of δE nadd (SIBFA) as compared to the situation at the HF level. A larger magnitude for ∆E(DFT) value of -79.1 kcal/mol was recently reported by Paragi and Fonseca Guerra 36 but this was possibly due to the use of different basis set and functional (BLYP-D/TZ2P) and/or stacked geometry. b). Stacked half-quartets. The results are reported in Table 3 and 4 at HF/RVS and correlated levels, respectively, regarding the four-guanine complex and the six separate dimeric complexes G1-G2 till G1'-G2'.

Correlated results (Table

HF/RVS calculations. ∆E(HF) is overestimated by ∆E(SIBFA) by 4 kcal/mol out of 8. This was unexpected, the more so, as it is due to E MTP* accountable for a 5.8 kcal/mol overestimation, only partly compensated for by the other contributions. Examination of the individual dimers shows the error in E MTP* to originate mainly from the three stacked complexes G1-G2, G2-G1', and G2-G2', namely 1.5 kcal/mol for each. This could again diagnose some shortcoming in the present derivation of the distributed multipoles, more acute in stacked arrangements where the quadrupolar contributions could be more effective than in the H-bonded complexes. Regarding the second-order contributions, we observe δE nadd to now be positive, indicative of anticooperativity. It is slightly more pronounced in SIBFA than in

ACS Paragon Plus Environment

The Journal of Physical Chemistry QC/RVS (2.2 versus 1.7 kcal/mol). It stems predominantly from E pol , nevertheless E ct has a small, yet beyond rounding-off error, contribution to it.

Correlated calculations. ∆E tot (SIBFA) of -38.0 kcal/mol now compares much more favorably to ∆E(B3LYP-D3) and ∆E(B97-D3) of -39.8 and-34.8 kcal/mol, respectively. The agreement is better with the former functional, a reverse situation to the one occurring with the H-bonded complexes. There is a significantly more uniform agreement with the six guanine dimers than was the case at the HF level. This could indicate a better balanced multipolar distribution derived at the DFT than at the HF level. But it is not precluded that the choice of the parameters used in the GDMA procedure was better adapted for the DFT-derived multipoles than those used for the HF ones. δE nadd of 1.5 kcal/mol for SIBFA is closer to δE nadd for B3LYP-D3 (1.6 kcal/mol) than for B97-D3 (2.2 kcal/mol).

In the absence of metal cations, such values of δE nadd are about five times smaller than those found with an opposite sign for each of the two quartets, so that the stacking of two G quartets should be dominated by cooperativity. However it should be noted that there are four equivalent 'half-quartets', G1-G2/G1'-G2' till G4-G1/G4'G1', each one sharing two G bases from the previous one. A complete analysis would need to be done involving all eight bases of the two stacked GQ complexes, with as well as without metal cations. It is planned in future studies. At this stage it was essential to evaluate how well the trends and magnitudes of δE nadd found in QC could be retrieved in SIBFA. 5). The results are fully consistent with those found at the quartet level. At the HF level), ∆E(SIBFA) is larger than ∆E(HF) by an amount nearly twice as large as with one quartet (the BSSE correction was not subtracted). At the correlated level, ∆E tot (SIBFA) of -169.8 kcal/mol is larger in magnitude than ∆E(B97-D3) of -168. thus results into additional stabilization energies, δ stack , of -36.6, -40.8, and -30.3 kcal/mol in these respective approaches. The close δ stack values found in SIBFA and B97-D3 imply that the additional stabilization brought by stacking and its impact on cooperativity, are computed consistently by these two approaches. It is noteworthy that their difference, 4.2 kcal/mol, is smaller than the B97-D3 and B3LYP-D3 δ stack difference of 10.5 kcal/mol.

c). Two stacked G quartets (Table

We deem these preliminary analyses to be essential prior to studying the cation complexes of the stacked G quartets. With two cations thus ten interacting entities, there is a total of fortyfive 'dimer' interactions, namely twenty-eight base-base, sixteen base-cation, and one cationcation interaction. How well the H-bonded and stacked base-base interactions, modulated by the potentials and fields of the other bases and the cations, are accounted for by APMM methods, can critically impact the potential energy surfaces they generate for the stacked G quartet-cation complexes 2. Channeling of one monovalent cation (Li + , Na + , K + , Rb + ) through two stacked G quartets.

The cation is displaced along the Z axis from bottom to top in fifteen steps of 0.25 Å amplitude. We retain the coordinates adopted by Gkionis et al. [START_REF] Gkionis | Ion Binding to Quadruplex DNA Stems. Comparison of MM and QM Descriptions Reveals Sizable Polarization Effects Not Included in Contemporary Simulations[END_REF] . However the atom ordering is not the same. The one used in the SIBFA library of fragments starts with H(9), N9, C8, etc.

The G8 complex was reconstructed by molecular graphics to superimpose with the one given by Gkionis et al. With this construction, the planar separation between the two quartets is 3.4 Å, but while each quartet is flat, there can be small deviations from the uniform Z values of -1.7 and +1.7 Å for the bottom and the top quartets. The geometries are given in Supp. Info.

The cation starts at Z=-2.72 Å, 1 Å below the center of the bottom quartet, and ends at Z=0.78 Å, 2.5 Å above that center and 0.9 Å below the center of the top quartet. The energy profiles for the channeling of each cation are reported in Figures 2a-d Both ∆E(HF) and ∆E(SIBFA) have a shallow maximum at Z=0 Å, when the cation is inbetween the two quartets, with values of -206.7 and -202.5 kcal/mol, respectively. In SIBFA, E 1 has a much larger magnitude than E 2 , with a ratio in the range 1.4-1.65. However, while it is E 1 that confers its shape to ∆E(SIBFA) in the first half of the trajectory, namely until Z=-1, it is the reverse in the second half, and the shallow maximum at Z=0 Å is clearly due to E 2 .

Na+ (Figure 2b). Both ∆E(HF) and ∆E(SIBFA) curves are very shallow in-between Z=-1.22 and Z=0.78 Å, where the energy variations are <1.5 kcal/mol out of 205, i.e. less than 1%.

The curves are actually shallower than the corresponding Li + curve, despite the larger size of the cation. Two very flat minima can still be discerned. The first is at Z=-0.97 Å (HF) and at Z=-1.22 Å (SIBFA) with ∆E(HF) and ∆E(SIBFA) values of -205.6 and -207.2 kcal/mol, respectively. The second minimum is at Z=0.78, with virtually identical HF and SIBFA values of -206.2 kcal/mol. As was the case with Li + , the shape of ∆E(SIBFA) is dictated by E 1 in the first half of the trajectory and by E 2 in the second half. K + (Figure 2c). The HF and SIBFA energy profiles dramatically differ from the Li + and Na + ones, and they do so on many counts. There now is one single, well-defined minimum at Z=0, namely when the cation is in the center of the cavity (∆E(HF)= -196.3 and ∆E(SIBFA)=-195.8 kcal/mol). This is a clear indication for the fact that owing to its size, K + is the most suited cation to bridge the two G quartets at equal distances. Along the alkali cation series, it is also the cation endowed with the most favorable stacked G quartet binding affinity as recently demonstrated by QC energy balance analyses [START_REF] Zaccaria | The Role of Alkali Metal Cations in the Stabilization of Guanine Quadruplexes: why K+ is the Best[END_REF] . The relative weight of E 1 to E 2 fluctuates much more than with the two smaller cations: from 1.26 in the first third of the curve to 1.74 at, and close to the minimum.

Rb + (Figure 2d). The HF and SIBFA curves have much more pronounced features than the K + curves. Both have a well-defined maximum at Z=-1.72 Å, namely when the cation in the center of the first monomer, a clear indication of the fact that the cation is oversized with Page 11 of 55

ACS Paragon Plus Environment

The Journal of Physical Chemistry respect to the GQ cavity size. As for K + , the energy minima are at Z=0 Å, in the center of the cavity, but are deeper than for K + . ∆E(HF) and ∆E(SIBFA) at the minimum are numerically close (-186.7 and -190.0 kcal/mol, respectively), and this is the case for the entire trajectory, except to a limited extent (up to 4%) in the high-lying zone. The ratios of E 1 to E 2 are now very contrasting. E 1 dictates its shape to the entire trajectory. Except for the first point, it is lesser in magnitude than E 2 until Z=-1.2 Å, and has the least value at Z=-1.72 Å, when the cation is in the center of the first quartet.

b) ∆E(DFT-D3) and ∆E tot (SIBFA) profiles.

We will mostly display the results found with the B97-D3 functional, which in keeping with previous results was the one giving rise to the closest agreements with ∆E tot (SIBFA) [START_REF] El Hage | Substituent-Modulated Affinities of Halobenzene Derivatives to the HIV-1 Integrase Recognition Site. Analyses of the Interaction Energies by Parallel Quantum Chemical and Polarizable Molecular Mechanics[END_REF][START_REF] Gresh | Complexes of a Zn-Metalloenzyme Binding Site with Hydroxamate-Containing Ligands. A Case for Detailed Benchmarkings of Polarizable Molecular Mechanics/Dynamics Potentials when the Experimental Binding Structure is Unknown[END_REF] . The ∆E(B3LYP-D3) curves, shown below for the Li + complexes, are invariably parallel to the B97-D3 one but app. 30 kcal/mol lower. The increases in magnitude upon passing from ∆E(HF) to ∆E(B97-D3) and from ∆E(SIBFA) to ∆E tot (SIBFA), are closely similar for all four cations. E disp (SIBFA) stems solely from the guanine-guanine interactions contributing -88.5 kcal/mol for all complexes. It is not shown in the figures.

Li + (Figure 3a). While ∆E tot (SIBFA) has the same shallow profile as ∆E(SIBFA) with a flat minimum at Z=-1.22 Å. ∆E(B97-D3) undergoes a regular decrease until Z= -0.47 Å. The decrease then becomes more accented with a minimum at Z=0.03 Å, and this is observed as well for ∆E(B3LYP-D3) but had not happened with ∆E(HF). This appears to stem from the D3 functional. Thus Figure 3b recasts the profile computed with the B3LYP functional in the absence of such correction, seen to parallel ∆E(HF), but not ∆E(B3LYP-D3).

Na + (Figure 3c). ∆E(B97-D3) has a profile similar to the Li + one, with its minimum at Z=0.03, while ∆E tot (SIBFA) has a shallow behavior from Z=-1.7 Å throughout, with two shallow minima at Z=-1.22 Å and Z=0.78 Å. As a result, the relative error of ∆E tot (SIBFA) with respect to ∆E(B97-D3) is the smallest at Z=-1.72 Å, at the center of the first quartet (2%), and is the largest at Z= 0 Å (6%). K + (Figure 3d). ∆E(B97-D3) and ∆E tot (SIBFA) have similar shapes, as in the absence of correlation/dispersion (Figure 2c). There is a shallow increase in magnitude in the region Z= -2.72 till -0.97 Å, followed by a more accented descent to the minimum at Z=0 Å. As with Na + , the minimum is deeper for ∆E(B97-D3) than for ∆E tot (SIBFA), and the relative error is the largest there (4%) while it was <1% at the HF level.

Rb + (Figure 3e). ∆E(B97-D3) and ∆E tot (SIBFA) have more parallel shapes than was the case with the smaller cations, the offset never exceeding 11.3 kcal/mol ouf of 220, its average value being of 7.3 kcal/mol, the average energy difference amounting to 3%. As with the HF calculation, the maximum of both ∆E(B97-D3) and ∆E tot (SIBFA) is at -1.72 Å in the center of the first quartet, and the minimum is at Z=0.0. We report in Appendices III and IV the Tables giving the energies plotted in Figures 4a-d and5a-d. We do not list the separate values of E pol and E ct within E 2 . E ct has shallower variations than E pol , and is in the range -18 ---24 kcal/mol. a) ∆E(HF) and ∆E(SIBFA) profiles.

Channeling of two monovalent cations (Li

The results are reported in Figures 4a-d. For all cations, both QC and SIBFA profiles very significantly differ from the corresponding one-cation profiles.

Li+ (Figure 4a). The minima of both QC and SIBFA are shifted from Z=-1.47 to Z=-1.72 Å, namely in the center of the first quartet, after which the energy raises continuously rather than plateauing. The E 1 /E 2 ratio decreases progressively, until E 2 becomes larger in magnitude at Z=-0.22 Å. This clearly translates the raise in total polarization energy of the two G quartets upon shortening the cation-cation distance, counteracting the raise in their mutual Coulomb repulsion. It is notable that the sum of E 1 and E 2 enables SIBFA to match ∆E(HF) with relative errors < 2% thus accounting for the Z-dependent ∆E magnitude increase upon passing from the one-to the two-cation G8 complexes.

Na + (Figure 4b). The profiles of ∆E(HF) and ∆E(SIBFA) both resemble the Li + ones, with a minimum shifted at Z=-1.72 Å, a continuous raise thereafter, and a E 1 /E 2 crossing at Z=-0.22.

In the low-energy zone, ∆E(SIBFA) matches ∆E(HF) with relative errors of 2%, comparable to those found with Li + , but now overestimates it rather than underestimating it. K + (Figure 4c). The Z=0 Å minimum has disappeared owing to the increased electrostatic repulsion undergone by the moving cation. Instead, the energy plateaus along the whole Z=-2.72 ---0.22 Å region, undergoing <2% variations. Except for the first two points, E 2 is now larger in magnitude than E 1 , the gap increasing steeply past Z=0 Å. The relative error of ∆E(SIBFA) with respect to ∆E(HF) is 3%.

Rb + (Figure 4d). With respect to the K + curve, a shallow maximum is present at Z=-1.72 Å, when the moving cation is in the center of the first quartet, and there still does exist a minimum at Z=-0.47 Å, 0.5 Å before the center of the cavity is reached. Past this minimum, the energy raises more steeply than with K + , owing to the much larger short-range Rb + -Rb + than K + -K + short-range repulsion. Thus with both QC and SIBFA there are two well-defined minima, at the start of the trajectory (Z=-2.72 Å) and close to the center of the cavity (Z=-0.47 Å). E 2 is now invariably larger in magnitude than E 1 , but the shape of ∆E is dictated by E 1 .

The relative errors of ∆E(SIBFA) are larger than in the two K + case, and can reach 6% in the low-energy zone.

b) ∆E(B97-D3) and ∆E tot (SIBFA) profiles.

Li+ (Figure 5a). Both QC and SIBFA curves closely resemble their counterparts at the HF level (Figure 4a). There is a shallow minimum at Z=-1.72 Å, in the center of the first quartet followed by a progressive rise in the energy. The E 1 and E 2 curves cross at Z=-0.72 Å, 0.30 Å earlier than at the HF level. The relative SIBFA versus QC error now averages 5%, larger than the 2% error found at the HF level.

Na + (Figure 5b). Again, the shapes of ∆E(B97-D3) and ∆E tot (SIBFA) are very similar to their HF counterparts. The shallow minimum in the first half of the trajectory locates Na + at Z=-1.47 Å, namely 0.25 Å above the center of the first quartet, and the E 1 /E 2 crossing occurs 0.5 Å earlier than at the HF level. The relative SIBFA error averages 3% over the relevant energy zone but is lowest at the minimum. K + (Figure 5c). Both curves plateau closely similar to their HF counterparts. There is a discernible minimum for Z=-0.72 Å, and the total energies start to rise steeply as the moving cation progresses 0.5 Å past it and beyond. The crossing of the E 1 /E 2 curves now occurs at the very start of the trajectory, rather than 0.5 Å past it as with the HF curves. The relative errors in SIBFA are strikingly small (1%) over the low energy range of distances, namely in thirteen points covering a 3 Å trajectory.

Rb + (Figure 5d). As with the HF curves, a local maximum is located in the center of the first quartet, and for both SIBFA and QC the global minimum occurs at Z =--0.47 Å, 0.5 Å before the center of the cavity. E 2 has larger magnitudes than E 1 throughout and the gap increases along the entire trajectory except close to the energy minimum. The relative errors in SIBFA are <2% for the thirteen first points in the trajectory.

Non-Covalent Interaction (NCI) analyses.

Analyzing the sign of the second eigenvalues, λ 2 , of the electron density Hessian matrix shows the nature of the non-covalent interaction: attractive if λ 2 is negative, and repulsive if positive. If ρ denotes the electronic density, a plot of the reduced gradient density, s, as a function of the product λ 2 ρ provides information on such a nature, classified as attractive, weak, or repulsive.

A red-green-blue color scheme is chosen for the 2D and 3D plots to distinguish between strong and attractive interactions such as ionic or hydrogen bonds (in blue), weak interactions such as van der Waals or related interactions (in green), and repulsive interactions (in red).

We first consider the two stacked quartets prior to cation binding (Figure 6). The onset of stabilizing inter-quartet interactions is attested by the large, deep-green isosurfaces in the midplane (Figures 6a-b). Within each quartet plane, the spherical blue-green isosurfaces between H(N1) and O6, and between H(N2) and N7, attest to the onset of strong inter-molecular Hbonds. The red oval surfaces in the middle of both rings of each guanine attest to local intramolecular repulsions (Figure 6a).

We next consider the NCI plots for one K + cation channeling along the Z axis (Figures 7a-d).

When it is located beneath the first quartet (Z=-2.7 Å, Figure 7a), the K + --O6 interactions are translated by downward-pointing green disk-shaped isosurfaces, with no visible perturbation of the mid-plane isosurfaces. Its in-plane location (Z=-1.7 Å, Figure 7b) translates by in-plane alignments of O6, K + , and the disk-shaped isosurfaces, with little impact on the mid-plane isosurfaces. The location of K+ in the center of the cavity (Z=0.0 Å, Figure 7c) translates, on the one hand, by an upward displacement of the disk isosurfaces connecting K + to the bottom quartet O6 and downward displacement of the disk isosurfaces connecting it to the top quartet O6. The mid-plane isosurfaces are less disjoint than with Z=-2.7 and -1.7 Å, or in the absence of cation altogether. When K + is closer to the top quartet (Z=0.7 Å, Figure 7d), the O6-K + isosurfaces stemming from the bottom quartet have faded out, while those stemming from the top quarter have their density increased. The mid-plane isosurfaces are not significantly impacted.

The NCI features along the four steps of this trajectory become more accented in the presence of an additional cation held on top of the top quartet at Z=2.8 Å (Figures 8a-d). Whatever the location of the channeling cation, there are four upward directed disk-shaped isosurfaces between the four O6 atoms of the top quartet and the topmost cation. As was the case with one K + , the corresponding four isosurfaces of the bottom quartet follow the moving cation: downward-pointing for Z=-2.8 Å (8a), in-plane location for Z=-1.8 Å (8b), upward-pointing for Z=0.0 Å (8c), and fading out for Z=0.8 Å (8d). For Z=0.0 Å, it is instructive to observe that each O6 atom from the top quartet contributes two disk-shaped isosurfaces: upward-pointing toward the top K + , and downward-pointing toward the central K + . The mid-plane isosurfaces are brighter than the corresponding ones with only one K + .

The different peaks present in the 2D plot can be analyzed for a better understanding of the features of the 3D isosurfaces. We can detect first the signature peak of the guanines, prior to K + binding (Figure 6c). There is a neat peak in the attractive region (λ 2 ρ =-0.02 au). It corresponds to the inter-molecular H1-O6 and H2-N7 bonds. The two green peaks (λ 2 ρ = -0.005 and 0.0 au) correspond to the mid-plane isosurfaces. The red peak at λ 2 ρ=0.025 au corresponds to the repulsive intra-ring interactions in each guanine.

We next consider the evolutions of these peaks in the presence of one and two K + cations.

When K + is below the bottom quartet (Z=-2.8 Å, Figure 7a), an additional blue peak appears at λ 2 ρ =-0.015 au. When K + is in the center of this quartet (Z=-1.8 Å, Figure 8b), the peak is shifted to a more attractive region with λ 2 ρ =-0.03 au. When K + moves to the center of the cavity (Z=0.0 Å, Figure 7c), or closer to the topmost quartet (Z=-0.8 A, Figure 7d), this blue peak disappears. For the first two points (Z=-2.8 and -1.8 Å) and for the last one (Z=+0.8 Å, Figure 7d), the two central green peaks on both sides of λ 2 ρ =0 have the same features as in the absence of K + . When K + is in the center of the cavity (Z=0.0 Å), the left-most peak becomes wider, consistent with the increase of the mid-point isosurfaces in the 3D plots (Figure 7c).

We next consider the two K + complexes. There are no marked differences for the first two points (Z=-2.8 and -1.8 Å, Figures 8a-b) with respect to the one-cation case. When K + is in the center of the cavity, on the other hand, both left and right green peaks around λ 2 ρ= 0.0 au are enlarged with respect to the one-cation case. At Z=0.8 Å, an additional red peak appears at λ 2 ρ =0.015 au. It translates the onset of short-range repulsive interactions between the two cations which are at 2.0 Å from one another.

For both Na + and Rb + complexes, both 2D and 3D plots present very similar features as with K + (not shown). The most conspicuous feature which differentiates between them is found with the two-cation complexes at Z=0.7 Å. It relates to the repulsive peak at λ 2 ρ = 0.015 au, more spread-out with Rb + than with K + , and virtually undetectable with Na + .

In conclusion, the NCI analysis demonstrates that both the guanine-guanine and guaninemetal key interactions are in a van der Waals regime. The strength of guanine-metal interactions clearly depends on the cation position, being maximized when the metal is located at the center of the cavity. The relative selectivity of the various cations therefore Page 16 of 55

ACS Paragon Plus Environment

The Journal of Physical Chemistry relies on subtle variations of van der Waals interactions that can be highlighted by the presented supermolecular interaction energies.

Discussion.

We have carried out benchmark comparisons of the SIBFA APMM procedure against QC results in a demanding test case, the complex of two stacked guanine quartets (GQ) with one or two alkali cations Li + , Na + , K + and Rb + , and the energy profile for the channeling of one cation along the Z axis with or without another one set above the topmost quartet. These tests are motivated by the emerging importance of GQs in molecular biology, in drug design, and in nanoscience. The alkali cations K + and Na + inside their cavity can play a major role in their stabilization. Stacked GQ's in a wealth of settings clearly constitute a privileged avenue for a manifold of future large-scale APMM applications. They are stabilized by the simultaneous interplay of multiply H-bonded interactions, stacking interactions between extended and polar surfaces, and cation-ligand interactions. Non-additivity effects strongly come into play, and have to be correctly quantified by APMM procedures. It is thus mandatory to objectively benchmark their expectable accuracy against high-level QC, and how well they perform compared to non-polarizable potentials [START_REF] Gkionis | Ion Binding to Quadruplex DNA Stems. Comparison of MM and QM Descriptions Reveals Sizable Polarization Effects Not Included in Contemporary Simulations[END_REF] . We have adopted a staged approach comparing first ∆E(SIBFA) to ∆E(HF) and then ∆E tot (SIBFA) to ∆E(DFT-D3). An essential motivation for this is the fact that all ∆E(SIBFA) parameters were derived beforehand on the basis of RVS EDA's done at the HF level with the aug-cc-pVTZ(-f) level [START_REF] Devereux | A Supervised Fitting Approach to Force Field Parametrization with Application to the SIBFA Polarizable Force Field[END_REF] , then extended to alkali cations [START_REF] Dudev | Quantum-Chemistry Based Calibration of the Alkali Metal Cation Series (Li+-Cs+) for Large-Scale Polarizable Molecular Mechanics/Dynamics Simulations[END_REF] and to cytosine and guanine [START_REF] Gresh | Analysis of the Intermolecular Interaction Energies by Parallel Quantum Chemistry and Polarizable Molecular Mechanics[END_REF] . The passage to correlated tests is done solely upon replacing the HF-derived multipoles and polarizabilities by those derived at the DFT level, and adding an appropriately rescaled E disp contribution. Can we, then, be equally successful, if at all, in reproducing in succession both HF and DFT-D calculations?

We considered the stabilization energies of one, then of two stacked quartets prior to cation binding, and then the complexes of two stacked GQ's with one and with two cations.

-In-plane G quartet. At the HF level, in each of the four doubly H-bonded dimers G1-G2 till G4-G1, E MTP and E rep were each accountable for overestimations of 0.5 kcal/mol out of 16 (E C ) and 11 (E exch ). This indicates that improvements to E MTP could be sought for, such as GEM or ISA. These could impact E rep, as well, since E rep embodies a dependence upon electrostatics in its prefactor. The second-order contributions, E pol and E ct gave close matches to their QC counterparts. Cooperativity was correctly accounted for, δE nadd amounting to -9.6 and to -9.2 kcal/mol in HF and SIBFA, respectively. At the correlated level, ∆E tot (SIBFA) of -68.8 kcal/mol was found to have a closer agreement to ∆E(B97-D3) of -69.8 kcal/mole than to ∆E(B3LYP-D3) of -74.7.

-Two stacked half-quartets. Such an arrangement involves a doubly H-bonded dimer of the first quartet and one of the two most overlapping H-bonded dimers of the other quartet. δE nadd now has a positive sign, indicative of local anticooperativity, and its QC value is satisfactorily accounted for in SIBFA: 1.5 vs. 2.2 kcal/mol at the HF level, and 1.5 versus 1.6 and 2.2 kcal/mol at the B3LYP-D3 and 2.2 levels, respectively. As with the in-plane G quartet, better agreements are found between ∆E tot (SIBFA) and ∆E(DFT-D3) than at the HF level. Thus ∆E tot (SIBFA) amounts to -38.0 kcal/mol as compared to -39.8 and -34.8 kcal/mol with the B3LYP-D3 and B97-D3 functionals.

-Two-stacked GQ complex. The best agreement was between ∆E tot (SIBFA) and ∆E(B97-D3), namely -169.8 vs. -168.6 kcal/mol. This value is smaller than the -184 kcal/mol ∆E(B3LYP-D3) one. The agreement between SIBFA and B97-D3 is again closer than that of the two DFT-D3 procedures between themselves. This was observed previously in a Zn(II)-metalloprotein recognition site [68], and is also found in all mono-and dication curves.

-Cation channeling. We next compared the ∆E(SIBFA) and ∆E(HF) profiles for one cation channeling in the absence of the other. SIBFA could closely match the QC profiles for all four cations, the relative ∆E error being always <2% in the low-energy zones. Li + and Na + had similar, very flat profile, nevertheless two minima were discernible for both: the first, about 0.50 Å above the center of the first quartet, and the second, at the end of the trajectory. A flat maximum was at the center of the cavity, found by both QC and SIBFA. In SIBFA, the magnitude of ∆E was mostly due to E 1 accounting for its two thirds, yet it was E 2 that conferred to ∆E its shape in the last part of the trajectory. K + had a starkly different profile, increasing regularly in magnitude until it reached the center of the cavity. The E 1 /E 2 ratio was much more strongly variable than with the two cations, from 1.25 till 1.80 in the center of the cavity. The contrasted E 1 and E 2 profiles from Z=-1.47 Å and beyond are noteworthy, yet they sum up to confer the right profile to ∆E(SIBFA) as compared to ∆E(HF). A significant reshape of the energy profile, imposed by E 1 , occurred with Rb + , with an accented maximum in the center of the first quartet, and a minimum in the center of the cavity.

There are lesser agreements upon comparing ∆E tot (SIBFA) and ∆E(B97-D3) for the two lighter cations, Li + and Na + . ∆E(B97-D3) and ∆E(B3LYP-D3) contrary to ∆E(HF) now have a minimum at Z=0 Å. It is not found for ∆E tot (SIBFA), which behaves exactly as ∆E(HF) and ∆E(SIBFA), but also as ∆E(B3LYP). The reasons for the different B3LYP and B3LYP-D3 shapes remain to be clarified. ∆E tot (SIBFA) compares much more favorably to ∆E(B97-D3) with K + , the Z=0 Å minimum at the center of the cavity being the 'natural' one. Nevertheless the relative SIBFA error is still the largest at this point: 4% as compared to 2% in the first half of the trajectory. There is a more regular SIBFA/QC parallelism with Rb + , the relative error being contained to <3% throughout.

The energy profiles in the presence of a second cation are markedly different. Could, upon decreasing the cation-cation distance, the increase of magnitude of E 2 due to increased polarization correctly compensate for the decrease of E 1 magnitude due to increased cationcation electrostatic repulsion, and also, for the larger cations close to the end of their trajectory, the increase of their short-range repulsion? Any E 2 /E 1 imbalance, that could be further magnified by the non-additivities of E pol and E ct in the presence of the second cation, could severely downgrade the SIBFA versus QC parallelism. Regarding ∆E(SIBFA) and ∆E(HF), the lighter two cations had closely similar energy profiles, their minima now occurring in the center of the first quartet, rather than above it. The relative error is for both <2% over the whole low-energy zone. A cross-over of the E 1 and E 2 curves occurs at 0.25 Å before reaching the center of the cavity. A close parallelism is also found for K + . With both QC and SIBFA, the energy minimum is considerably shallower than in the one-cation case, and is displaced 0.75 Å beneath the center of the cavity. The E 1 /E 2 crossover now occurs already at 0.5 Å after the start of the trajectory. The relative SIBFA/QC error is < 3% in the low-energy range. The shape of the K + energy profile is further accented with Rb + . While the first maximum is still in the center of the first quartet, the minimum is displaced 0.5 Å beneath the center of the cavity. E 2 is larger in magnitude throughout the entire trajectory, the gap increasing upon Rb + progression along the Z axis. ∆E(SIBFA) very closely reproduces ∆E(HF) in the first part of the trajectory, the relative error being <2%, increasing slowly to 5% at the energy minimum. Closely similar conclusions can be reached regarding ∆E(B97-D3) and ∆E tot (SIBFA). For both Li + and Na + , the two curves display good parallelism, with relative SIBFA errors < 4% (Li + ) and <3% (Na + ). The energy minima are in the center of the first quartet. The E 1 /E 2 crossover occurs for ∆E tot (SIBFA) 0.50 Å earlier than for ∆E(SIBFA), translating the larger relative weight of E 2 with respect to E 1 with correlated multipoles and polarizabilities. For K + and Rb + , there is a rather remarkable overlap of the SIBFA and QC trajectories throughout. The relative errors are now 1%. Nevertheless the fact that we recover a better agreement with two than with one cation implies there are some compensations of errors that we will seek to unravel in future work. The E 1 /E 2 crossover now occurs at the very start of the trajectories. For K + , the energy minimum is at Z=-0.72 Å, consistent with its ∆E(HF) location. For Rb + , with both ∆E(B97-D3) and ∆E tot (SIBFA), it is shifted 0.25 Å upward compared to K + , and is more accented.

The present results could be used to benchmark other polarizable potentials [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF][START_REF] Zhang | Polarizable Multipole-Based Force Field for Aromatic Molecules and Nucleobases[END_REF][START_REF] Lemkul | Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model[END_REF][START_REF] Harder | Atomic Level Anisotropy in the Electrostatic Modeling of Lone Pairs for a Polarizable Force Field Based on the Classical Drude Oscillator[END_REF] or semiempirical QC methods. For this purpose we give as Supp. Info S1-S4 the G09 input data used for the G8 complexes with two cations. The moving cation is given its coordinates at the start of the trajectory (Z=-2.72 Å), the Z coordinate of the fixed cation being set at 2.80 Å.

Concluding remarks.

The present SIBFA benchmarks appear convincing and should lend credence to prospective simulations of stacked G quartets and their metal cation complexes. Long-duration molecular dynamics (MD) on periodic GQ stacks and their complexes with metal cations are thus envisaged. They will benefit from the recent integration of the SIBFA potential and its gradients into the massively parallel Tinker-HP package 73 . This should enable to study their ion transport properties, the impact of molecular environment in different phases, such as liquid, gel, solid, or close to a conducting surface, etc. MD could also provide poses for postprocessing by QC to calculate their most salient electronic properties. GQs could be connected together by variable linkers, to yield supramolecular complexes of tunable properties. We are also pursuing the refinements of the SIBFA potential. On the one hand, alternatives to GDMA multipoles are being tested. On the other hand, in line with our previous work [START_REF] Gresh | Analysis of the Intermolecular Interaction Energies by Parallel Quantum Chemistry and Polarizable Molecular Mechanics[END_REF] , additional SAPT calculations [START_REF] Jansen | Comment on "Using Kohn-Sham Orbitals in Symmetry-Adapted Perturbation Theory To Investigate Intermolecular Interactions[END_REF][START_REF] Jeziorski | Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes[END_REF][START_REF] Misquitta | Dispersion Energy from Density-Functional Theory Description of Monomers[END_REF] on the complexes of nucleic acid bases with cation and water probes are being performed. This should enable a finer-tuning, now at the correlated level, of each of the five SIBFA contributions against its QC counterpart. There are simultaneous ongoing SIBFA refinements of the sugar-phosphate backbone of DNA and RNA which should pave the way for future APMM simulations of quadruplex DNA and their ligand complexes. These will be reported subsequently. G09 input data used for the G8 complexes with two cations (Li + , Na + , K + , or Rb + ). The moving cation is given its coordinates at the start of the trajectory (Z=-2.72 Å) the Z coordinate of the fixed cation being set at 2.80 Å.

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Figures captions:

Page

1 .

 1 Stabilization energies in the G quartet, the stacked half-quartet, and in two stacked G quartets. a) G quartet. Tables 1 report the results of the QC and SIBFA intermolecular interaction energies at the HF and correlated levels respectively.

  6 kcal/mol by 1 kcal/mol out of 169, while it was smaller than it by 1 kcal/mol out of 69 in one quartet. It is smaller in magnitude than ∆E(B3LYP-D3) of -184.0 by 14 kcal/mol out of 180, an amount only slightly larger than the double of the 6 kcal/mol out of 70 found with one quartet. It is noted that the SIBFA and B97-D3 values of -169 kcal/mol are now close to the -172.4 reported by Paragi and Fonseca Guerra for the stacked GQ dimer 36 . It is instructive to compute the amount of stabilization brought by stacking, upon comparing these ∆E tot values to twice the values found in the isolated quartets. Such values amount to -137.6, -139.6, and -149.4 in SIBFA, B97-D3 and B3LYP-D3, respectively. The stacking of the two G quartets ACS Paragon Plus Environment The Journal of Physical Chemistry

  regarding ∆E(HF) and ∆E(SIBFA) computed with HF-derived multipoles and polarizabilities, along with the E 1 and E 2 contributions of ∆E(SIBFA). Figures 3a-e report the corresponding profiles for ∆E(B97-D3) and ∆E tot (SIBFA) now with correlated multipoles and polarizabilities and the dispersion contribution. The fifth point in the plot corresponds to the cation in the plane of the bottom quartet, and twelfth point to the cation at mid-distance from the two quartets. We report in Appendices I and II the Tables giving the energies plotted in Figures 2a-d and 3a-e. We do not list the separate values of E pol and E ct within E 2 . E ct has shallower variations than E pol , and is in the range -15 ---20 kcal/mol. a) ∆E(HF) and ∆E(SIBFA) profiles. The four cations display markedly different QC energy profiles, which can be closely matched by SIBFA. Li + (Figure 2a) ∆E(HF) has two shallow minima, at Z=-1.47 Å, 0.25 Å above the first quartet, and at Z=0.78 Å, with values of -211.3 and -210.4 kcal/mol, respectively. ∆E(SIBFA) has minima at these two points as well, with values of -209.1 and -207.0 kcal/mol, respectively.

  + , Na + , K + , Rb + ) through two stacked G quartets.

Appendices 1 - 4 .Appendix 1 .

 141 Channeling of one cation (Li + -Rb + ) along the Z axis of two stacked G4 tetramers. Comparison of ∆E(QC) and ∆E(SIBFA) values. Consistent with Figures 2-5, the Z values (first column) are multiplied by a factor of 100 for clarity. Channeling of one cation along the Z axis of two stacked G4 tetramers. Comparison of ∆E(HF) and ∆E(SIBFA) values (in kcal/mol

30 of 55 ACSFigure 1 .

 551 Figure 1. Complex of two stacked G quartets with two monovalent cations. Representation of the trajectory followed by the moving cation. The fixed cation is located on top of the upper quartet.

Figures 2 .Figures 3 .

 23 Figures 2. Complex of two stacked G quartets with one monovalent cation. ∆E(HF), ∆E(SIBFA), and E 1 /E 2 (SIBFA) profiles for cation channeling along the Z axis. 2a: Li + ; 2b: Na + ; 2c: K + ; 2d: Rb + . Figures 3. Complex of two stacked G quartets with one monovalent cation. 3a: ∆E(B97-D3), ∆E(B3LYP-D3), ∆E tot (SIBFA) and E 1 /E 2 (SIBFA) profiles for Li + channeling along the Z axis. 3b: ∆E(HF), ∆E(B3LYP) and ∆E(SIBFA) profiles for Li + channeling. 3c-3d: ∆E(B97-D3), ∆E(B3LYP-D3), ∆E tot (SIBFA) and E 1 /E 2 (SIBFA) profiles for cation channeling along the Z axis. 3c : Na + ; 3d : K + ; 3e: Rb + .

Figure 4 .

 4 Figure 4. Complex of two stacked G quartets with two monovalent cations with one moving cation. ∆E(HF), ∆E(SIBFA), and E 1 /E 2 (SIBFA) profiles for cation channeling along the Z axis. 4a: Li + ; 4b: Na + ; 4c: K + ; 4d: Rb + .

Figure 5 .

 5 Figure 5. Complex of two stacked G quartets with two monovalent cations with one moving cation. ∆E(B97-D3), ∆E tot (SIBFA), and E 1 /E 2 (SIBFA) profiles for cation channeling along the Z axis. 5a: Li + ; 5b: Na + ; 5c: K + ; 5d: Rb + . Figures 6a-b. 3D plot NCI analysis of the two guanine quartets without cation: (a) in the XZ plane; (b) in the XY plane; and (c) corresponding 2D plot. In the 3D plots, the NCI surfaces correspond to s = 0.50 isosurfaces. Figures 7a-d. 2D and 3D plot NCI analyses of the two guanine quartets with K + along the Z axis: (a): Z = -2.7 Å, (b): Z = -1.7 Å, (c): Z = 0.0 Å, (d): Z = 0.8 Å. In the 3D plots, the NCI surfaces correspond to s = 0.50 isosurfaces.Figures 8a-d. 2D and 3D plot NCI analyses of the two guanine quartets with a fixed K + cation at Z = 2.8 Å and a mobile K + cation along the Z axis : (a): Z = -2.7 Å, (b): Z = -1.7 Å, (c): Z= 0.0 Å, (d): Z= 0.8 Å. In the 3D plots, the NCI surfaces correspond to s = 0.50 isosurfaces.

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Representation of the complex of two stacked guanine quartets with one monovalent cation channeling along the Z axis and a fixed cation set on top of the upper quartet.

Figure 4 .Figure 5 . 51 Figure 6 .

 45516 Figure 4. Complex of two stacked G quartets with two monovalent cations with one moving cation. ∆E(HF), ∆E(SIBFA), and E 1 /E 2 (SIBFA) profiles for cation channeling along the Z axis. 4a: Li + ; 4b: Na + ; 4c: K + ; 4d: Rb + . For clarity, the Z values along the abscissa (Å) were multiplied by a factor 100.

  ).Channeling of one cation along the Z axis of two stacked G4 tetramers. Comparison of ∆E(∆E(B97-D3) and ∆E(SIBFA) values (in kcal/mol).Channeling of one cation along the Z axis of two stacked G4 tetramers in the presence of a fixed cation above the top tetramer. Comparison of ∆E(B97-D3) and ∆E tot (SIBFA) values (in kcal/mol).

		Li + Li +		Na + Na +		K + K +		Rb + Rb +
		∆E(QC) ∆E tot (SIBFA)	∆E(QC) ∆E tot (SIBFA)	∆E(QC) ∆E tot (SIBFA)	∆E(QC) ∆E tot (SIBFA
		∆E(QC) ∆E tot (SIBFA)	∆E(QC) ∆E tot (SIBFA)	∆E(QC) ∆E tot (SIBFA)	∆E(QC) ∆E tot (SIBFA
	-272	Li + -327.72	-309.20	Na + -317.50	-307.44	K + -288.21	-284.77	Rb + -259.11	-259.59
	-272 -247 -247 -222	∆E(HF) ∆E(SIBFA) -271.74 -256.79 -332.32 -315.79 -276.48 -264.35 -336.50 -321.28	∆E(HF) ∆E(SIBFA) -265.75 -255.82 -321.04 -313.30 -270.44 -262.68 -323.81 -318.00	∆E(HF) ∆E(SIBFA) -250.46 -243.78 -287.36 -285.81 ∆E(HF) ∆E(SIBFA) -232.92 -230.18 -253.71 -253.69 -250.82 -245.82 -231.08 -225.29 -285.65 -285.69 -247.81 -246.26
	-222 -197	-281.91 -339.33	-270.94 -325.19	-274.47 -325.72	-268.51 -321.26	-250.45 -284.13	-246.84 -285.24	-226.75 -243.16	-219 -239.96
	-197 -172	-286.17 -340.66	-276.10 -327.17	-277.8 -326.85	-273.07 -322.97	-250.48 -283.81	-247.68 -285.36	-223.98 -241.43	-213.99 -237.61
	-172 -147	-289.17 -340.58	-279.56 -327.11	-280.57 -327.32	-276.30 -323.09	-252.05 -285.12	-249.33 -286.53	-224.43 -243.54	-213.15 -240.58
	-147 -122	-291.05 -338.98	-281.25 -325.13	-283.02 -327.02	-278.23 -321.69	-255.67 -287.71	-252.29 -288.41	-229.10 -248.87	-217.9 -247.80
	-122 -97	-291.93 -335.82	-281.37 -321.58	-285.21 -325.79	-278.98 -318.89	-261.08 -290.62	-256.33 -290.09	-237.11 -255.40	-227.31 -256.49
	-97 -72	-292.05 -331.35	-280.35 -316.92	-287.12 -323.73	-278.78 -314.93	-267.37 -292.47	-260.65 -290.54	-247.40 -260.86	-238.73 -263.54
	-72 -47	-292.21 -326.14	-278.78 -311.64	-288.94 -320.85	-277.98 -310.14	-273.23 -292.27	-264.4 -289.03	-256.94 -262.92	-249.35 -266.54
	-47 -22	-295.15 -320.40	-277.29 -306.06	-290.76 -316.51	-277.07 -304.74	-277.7 -288.85	-267.16 -284.86	-264.26 -259.13	-257.39 -263.39
	-22 3	-300.27 -313.96	-276.41 -300.22	-292.14 -309.96	-276.53 -298.75	-280.09 -280.67	-268.87 -276.73	-268.93 -246.18	-262.23 -250.49
	3 28	-301.07 -306.29	-276.50 -293.73	-292.61 -300.67	-276.71 -291.74	-280.73 -265.66	-269.68 -261.19	-270.46 -216.46	-263.89 -218.52
	28 53	-299.39 -296.61	-277.65 -285.70	-292.07 -288.15	-277.70 -282.72	-279.85 -237.32	-269.64 -228.99	-268.25 -153.92	-262.42 -142.35
	53 78	-294.47 -283.64	-279.69 -274.66	-290.71 -269.87	-279.33 -269.97	-277.15 -178.80	-268.66 -154.97	-260.92 -21.58	-257.73 46.69
	78	-292.62	-282.20	-289.04	-281.20	-272.93	-266.59	-255.41	-249.91
		Supporting Information available for publication.

Table 1 . Intermolecular interaction energies (kcal/mol) and their contributions in a G4 tetramer and evaluation of non-additivities. Comparisons between

 1

	∆ ∆ ∆ ∆E(RVS) and

Page 32 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry

	Dimers	QC SIBFA
	G1-3, G2-4
	E c /E MTP	-1.3	-1.3
	E exch /E rep	0.0	0.0
	E 1	-1.3	-1.3
	E pol (RVS)/E pol *	-0.1	-0.1
	E pol (VR)/E pol	-0.1	-0.1
	E ct	0.0	0.0
	E ct *	0.0
	Summed values.
		-65.8 -67.8
	E c /E MTP
	Ee xch /E rep	44.4	42.4
	E 1	-21.4 -25.5
	E pol (RVS)/E pol *	-11.6 -11.4
	δE nadd	-5.9	-5.8
	E pol (VR)/Epol	-15	-13.4
	δE nadd	-8.8	-8.4
	E ct	-7.3	-8.8
	δE nadd	-0.8	-0.8

Page 33 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry

Table 2 . Intermolecular interaction energies (kcal/mol) and their contributions in a G4 tetramer and evaluation of non-additivities. Comparisons between

 2

	∆ ∆ ∆ ∆E(DFT-D3) and

Page 34 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry

Table 3 . Intermolecular interaction energies (kcal/mol) and their contributions in a stacked complex of bases G1-2 and G1'-2' of the first and second G4 tetramers and evaluation of non-additivities. Comparisons between

 3

	∆ ∆ ∆ ∆E(RVS) and ∆ ∆ ∆ ∆E(SIBFA).

Page 35 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry

	E ct	0.0	0.0
	E ct *	0.0
	∆E	-1.4	-1.3
	G1-G2'
	E c /E MTP	2.0	0.5
	E exch /E rep	3.0	3.6
	E 1	5.0	4.1
	E pol (RVS)/E pol *	-0.6	-0.7
	E pol (VR)/E pol	-0.6	-0.7
	E ct	-0.5	-0.2
	E ct *	-0.1
	∆E	3.9	3.2
	G2-G1'
	E c /E MTP	2.0	0.4
	E exch /E rep	3.1	3.7
	E 1	5.1	4.1
	E pol (RVS)/E pol *	-0.6	-0.7
	E pol (VR)/E pol	-0.6	-0.7
	E ct	-0.4	-0.2
	E ct *	-0.1
	∆E	4.0	3.2
	G2-G2'
	E c /E MTP	2.5	1.0
	E exch /E rep	2.3	2.9
	E 1	4.8	3.9
	E pol (RVS)/E pol *	-0.6	-0.6
	E pol (VR)/E pol	-0.6	-0.6
	E ct	-0.2	-0.1
	E ct *	-0.1

Page 36 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry

	∆E	3.8	3.2
	Summed dimer interactions	QC SIBFA
	E c /E MTP	-26.5 -32.0
	E exch /E rep	30.6	31.5
	E 1	4.1	-0.5
	E pol (RVS)/E pol *	-7.7	-7.7
	δE nadd	1.3	1.7
	E pol (VR)/E pol	-9.5	-8.7
	δE nadd	1.3	2.0
	E ct	-5.1	-4.8
	δE nadd	0.3	0.2
	E ct *	-3.1
	δE nadd	0.6	0.2
	∆E	-9.5	-14.0
	δE nadd	1.7	2.2

Page 37 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry

Table 5 . Intermolecular interaction energies (kcal/mol) in two stacked G4 tetramers.

 5

		QC(HF)	SIBFA
	∆E	-89.5	-94.8
		QC(B3LYP-D3)	QC(B97-D3)	∆E tot (SIBFA)
	∆E	-184.0	-168.6	-169.8

Page 39 of 55 ACS Paragon Plus Environment The Journal of Physical Chemistry 40

ACS Paragon Plus EnvironmentThe Journal of Physical Chemistry

Acknowledgments.

We wish to thank the Grand Equipement National de Calcul Intensif (GENCI), Institut du Développement et des Ressources en Informatique Scientifique (IDRIS), Centre Informatique de l'Enseignement Supérieur (CINES, France, project No. x2009-075009), and the Centre de Ressources Informatiques de Haute Normandie (CRIHAN, Rouen, France, project 1998053). JES and JS acknowledge support by the Czech Science Foundation [grant number 16-13721S]. Their research has also been financially supported by the Ministry of Education, Youth and Sports of the Czech Republic under the project CEITEC 2020 (LQ1601).

ACS