Cation competition and recruitment around the c-kit1 G-quadruplex using polarizable simulations.

Nucleic acid-ion interactions are fundamentally important to the physical, energetic, and conformational properties of DNA and RNA. These interactions help fold and stabilize highly ordered secondary and tertiary structures, such as G-quadruplexes (GQs), which are functionally relevant in telomeres, replication initiation sites, and promoter sequences. The c-kit proto-oncogene encodes for a receptor tyrosine kinase and is linked to gastrointestinal stromal tumors, mast cell disease, and leukemia. This gene contains three unique GQ-forming sequences that have proposed antagonistic effects on gene expression. The dominant GQ, denoted c-kit1, has been shown to decrease expression of c-kit transcripts, making the c-kit1 GQ a promising drug target. Toward disease intervention, more information is needed regarding its conformational dynamics and ion binding properties. Therefore, we performed molecular dynamics simulations of the c-kit1 GQ with K+, Na+, Li+, and mixed salt solutions using the Drude-2017 polarizable force field. We evaluated GQ structure, ion sampling, core energetics, ion dehydration and binding, and ion competition and found that each analysis supported the known GQ-ion specificity trend (K+ > Na+ > Li+). We also found that K+ ions coordinate in the tetrad core antiprismatically, whereas Na+ and Li+ align coplanar to guanine tetrads, partially because of their attraction to surrounding water. Further, we showed that K+ occupancy is higher around the c-kit1 GQ and its nucleobases than Na+ and Li+, which tend to interact with backbone and sugar moieties. Finally, we showed that K+ binding to the c-kit1 GQ is faster and more frequent than Na+ and Li+. Such descriptions of GQ-ion dynamics suggest the rate of dehydration as the dominant factor for preference of K+ by DNA GQs and provide insight into noncanonical nucleic acids for which little experimental data exist.

Significance

Here, we characterize the binding modes of different monovalent cations to the c-kit1 DNA G-quadruplex. We find that a balance of ion-G-quadruplex and ion-water interactions govern the kinetics of binding, coordination geometry, and ion competition around the nucleic acid structure. Our results help rationalize the preference for K+ ions over other monovalent cations and challenge the simplicity of the “optimal fit” hypothesis.

Introduction

Interactions between nucleic acids and ions underlie the physical, energetic, and conformational properties of DNA and RNA (1,2). Ion interactions influence nucleation of higher-order structure in nucleic acids as well as binding of proteins and drugs (1, 2, 3, 4). G-quadruplexes (GQs), highly stable noncanonical nucleic acid structures that function in transcription, translation, and telomere maintenance (5, 6, 7), are an important example of highly ordered structures governed by nucleic acid-ion interactions (8). Direct and solvent-mediated ion interactions modulate the folding of GQs and help maintain their stability (8, 9, 10). These stable structures are enriched in the promoter sequences of growth regulatory genes and proto-oncogenes and have been linked to neurodegenerative diseases (11), mental retardation (12), premature aging (13), and various types of cancer (14,15). Previous studies suggest that high-specificity ligands can be designed to modulate GQ stability and, ultimately, disease-associated gene expression (15, 16, 17, 18). Given their association with disease and dependence on monovalent cations (8), describing GQ-ion interactions is important for characterizing electrostatic interactions that are relevant to small-molecule drug design.

Topologically diverse GQs form in guanine-rich DNA and RNA, generally with sequences of GxNyGxNyGxNyGx (N is any nucleotide, x ≥ 3, and 1 ≤ y ≤ 7) (19). Interested readers are referred to (19) for a review of GQ structural features. These sequences self-associate into four-stranded helical structures composed of linker regions and a rigid core (20), which is formed from stacked tetrads of guanine bases (21). Each tetrad consists of a square planar arrangement of Hoogsteen hydrogen-bonded guanines, resulting in four inward-facing carbonyl oxygen atoms (O6). Because of this configuration, the central channel is strongly electronegative and requires cation stabilization (22). Monovalent metal ions are coordinated by tetrad guanines to counteract the repulsion arising from the proximity of the electronegative groups (22). This cation coordination is integral to GQ folding and is thought to contribute more to their stability than tetrad hydrogen bonding and stacking interactions (23).

The dominant cation found in GQ structures is K+, and this preference has implications for the kinetics and thermodynamics of GQ formation in vivo (24). Studies have found preferred folding pathways in the presence of K+ ions (25) and that in different buffer solutions, GQs can adopt different topologies, with variable thermal stabilities (26, 27, 28). Eisenman first reported GQ preference for K+ in 1962 and proposed that ion hydration was responsible for selectivity in GQs (29). However, he found that favorability of dehydration (Cs+ > Rb+ > K+ > Na+ > Li+) and melting temperatures of GQs followed different trends (K+ > Rb+ > Na+ > Cs+ > Li+). He also noted that K+ ions have an optimal fit to tetrad coordination sites (29,30). Based on these observations, Eisenman reasoned that selectivity for small cations is governed by a balance between electrostatic attraction and ion hydration as well as ion size (29,30). Since these pioneering studies, several experiments and techniques have supported Eisenman’s descriptions of GQ-ion binding. Using folding stability (31), thermal denaturation (32), relative ion binding (33), free energy of hydration (34), and van’t Hoff analysis (35), K+ preference over other cations has been confirmed consistently.

A quantum mechanical study of guanine tetrads and their ions further demonstrated K+ ions optimally fit in GQ tetrads. That is, K+ ions aligned symmetrically between stacked tetrads and were too large to adopt planar coordination, whereas Na+ and Li+ ions both adopted planar coordination with the tetrad (36). Although decades of studies emphasize the importance of ions in GQs, much is still unknown, including 1) the origin of K+ specificity, 2) how ions affect important loop regions and dictate folded topologies, and 3) the impact of ions on the electronic structure of the GQ core. These open questions reflect the unmet challenges of describing ionic atmospheres and their impact on the dynamic properties of nucleic acids (1,37). Though ions tightly bound to GQs can be observed by structural techniques (38), the ion atmosphere around GQs is highly dynamic and cannot be resolved using x-ray crystallography, NMR, or electron microscopy (1). Further, GQs are unstable in the absence of ions, making experimental analysis of ion binding difficult to perform. Considering these experimental limitations and that atomistic differences beyond the resolution of experimental studies may contribute to ion selectivity in GQs, molecular dynamics (MD) simulations are a promising approach to provide new insights.

Because GQ-ion selectivity is thought to be governed by a combination of electrostatic attraction, ion hydration, and charge density (29,30), MD simulations in explicit solvent are well suited to characterizing GQ-ion interactions. More specifically, polarizable force fields (FFs) should be employed for these simulations as the inclusion of explicit electronic polarization is necessary for accurately modeling bound ions in the GQ core (39). One such model for nucleic acids is the Drude-2017 polarizable FF (40,41), which attaches negatively charged auxiliary particles (Drude oscillators) to non-hydrogen atoms in the system via harmonic springs to model electronic degrees of freedom (42). These Drude oscillators reorient in response to changes in the local electric field that arise from electrostatic interactions within the system and geometric changes within each molecule (42). Details of the parametrization strategy, functional form, and integration algorithms associated with the Drude FF are reviewed in (42). The Drude FF has successfully modeled duplex DNA (40,41), DNA base flipping (43), ion-DNA interactions (44), DNA GQs (45, 46, 47), and RNA GQs (46), demonstrating its promise for investigations into the dynamics and energetics of GQs, including the c-kit1 GQ studied here.

The c-kit gene encodes for a receptor tyrosine kinase, and the folded c-kit1 GQ has been shown to decrease its expression (48,49). Aberrant c-kit gene expression is associated with gastrointestinal stromal tumors, mast cell disease, and leukemia (18,50, 51, 52, 53). As such, the structurally unique c-kit1 GQ (Fig. S1) is a promising target for high-specificity drug design (18,50,51). Previously, we studied the c-kit1 GQ using both the CHARMM36 (C36) additive (54) and Drude-2017 polarizable (40,41) FFs to contextualize polarization effects. In C36 simulations, we found issues with core-cation binding and preservation of loop structure, whereas the Drude FF provided improved descriptions of ion coordination in the tetrad core, reversible bulk ion binding, and structurally important loop interactions (45). Such improvements suggest that the inclusion of explicit electronic polarization is critical in studying GQs and that the Drude polarizable FF is a valuable tool to do so. Here, we employ the Drude-2017 FF to study the c-kit1 promoter GQ (53) in a variety of ionic conditions. These simulations build on our previous work and provide new insights into GQ interactions with K+, Na+, and Li+, competition among these monovalent cations, the ion atmosphere around GQs, and the influence of ion binding on the energetics and dynamics of GQs.

Methods

System construction

The starting c-kit1 GQ structure for all simulations was taken from the first model of the NMR ensemble in Protein Data Bank entry PDB: 2O3M (53). In the c-kit1 GQ, guanines 2, 6, 10, and 13 comprise tetrad 1; guanines 3, 7, 21, and 14 comprise tetrad 2; and guanines 4, 8, 22, and 15 comprise tetrad 3 (Fig. S1). Systems were constructed with bulk ions and coordinated ions in the tetrad core (“bound and bulk”) or with bulk ions and without coordinated ions (“bulk”). For “bound and bulk” systems, two ions were added between tetrads in bipyramidal antiprismatic coordination by assigning their positions as the average coordinates of the carbonyl oxygen (O6) atoms of the guanine bases of consecutive tetrads with the CHARMM program (55). The ion placed between tetrad 1 and tetrad 2 is referred to as “core ion 1,” and the ion placed between tetrad 2 and tetrad 3 is referred to as “core ion 2.” “Bulk” systems were constructed without bound ions to investigate bulk ion binding and the resulting dynamics of the c-kit1 GQ more completely.

All systems were initially constructed, minimized, and equilibrated using the C36 FF for nucleic acids (54,56,57) before being converted to the Drude-2017 nucleic acid FF (40,41), as described in our previous work (45). The c-kit1 GQ was centered in a cubic unit cell with a minimal box-solute distance of 10 Å, filled with the SWM4-NDP (58) water and a total salt concentration of ∼150 mM, including neutralizing counterions. The ion parameters of Yu et al. were applied, along with specific nonbonded parameters for nucleic acid-ion interactions (59,60). In total, 84 independent simulations were performed and are listed in Table 1.

Table 1

List of the contents and sizes of all simulation systems

System name	Solution	Starting conditions	Replicates	Box size (Å)
bound and bulk K+	∼150 mM KCl	two bound cations	3 × 1 μs	48
bound and bulk Na+	∼150 mM NaCl	two bound cations	3 × 1 μs	48
bound and bulk Li+	∼150 mM LiCl	two bound cations	3 × 1 μs	48
bulk K+	∼150 mM KCl	no bound cations	3 × 1 μs	48
bulk Na+	∼150 mM NaCl	no bound cations	3 × 1 μs	48
bulk Li+	∼150 mM LiCl	no bound cations	3 × 1 μs	48
extended mixed K+/Na+	∼75 mM KCl and ∼75 mM NaCl	no bound cations	3 × 1 μs	98
extended mixed K+/Li+	∼75 mM KCl and ∼75 mM LiCl	no bound cations	3 × 1 μs	98
mixed K+/Na+	∼75 mM KCl and ∼75 mM NaCl	no bound cations	40 × 200 ns	48
mixed K+/Li+	∼75 mM KCl and ∼75 mM LiCl	no bound cations	20 × 200 ns	48

Extended mixed systems

To characterize cation accumulation around GQs, a large box is needed, beyond what is required to avoid minimal image convention violations. Large systems allow for spatial resolution between ions in bulk solution and those that are accumulated around nucleic acids, thereby enabling a description of the ionic atmosphere around the GQ in terms of ion accumulation (61). These systems were prepared to allow us to describe the competition between different monovalent ions in terms of the excess ions that were accumulated around the c-kit1 GQ. Accordingly, the c-kit1 GQ was centered in a cubic box with an edge length of 98 Å, which was filled with SWM4-NDP water and equimolar mixtures of KCl and either NaCl or LiCl (Table 1).

Mixed systems

To assess the kinetics of ion binding, we set up multiple, short simulations in mixed salt solution. These simulations did not require large boxes as in the case of the extended simulations. Based on observations from our previous work with the c-kit1 GQ (45) and the extended mixed systems described above, we determined that 200 ns was sufficient time to observe a bulk ion binding event (defined as a bulk ion coming ≤3.5 Å from the center of tetrad 2). To perform statistical analysis on these observations, a large sample size was needed, so we prepared 40 replicates of mixed K+/Na+ and 20 replicates of mixed K+/Li+. The number of replicates for each system was determined by statistical power analysis, performed using R software (62).

MD simulations

Energy minimization, equilibration, and production runs were executed as described in our previous GQ simulation studies (45, 46, 47). Drude oscillator positions were relaxed with steepest descent minimization and adopted basis Newton-Raphson energy minimization in CHARMM. Then, NPT equilibration was carried out for 1 ns at 298 K and 1 atm by extended Lagrangian integration (63), implemented in NAMD as Langevin dynamics (63,64). Water and mobile ions were unrestrained during equilibration. The velocities of real atoms were scaled using an absolute thermostat at 298 K with a friction coefficient of 5 ps−1, and Drude oscillators were coupled to a low-temperature relative thermostat at 1 K with a friction coefficient of 20 ps−1. Periodic boundary conditions were applied in all spatial dimensions. The short-range van der Waals potential was switched to zero from 10 to 12 Å. All bonds to hydrogen atoms were constrained using the SHAKE (65) algorithm. The water molecules were kept rigid with SETTLE (66), allowing an integration time step of 1 fs. A “hard wall” constraint (67) was also applied to allow a maximal Drude-atom bond length of 0.2 Å to avoid polarization catastrophe. Equilibration was carried out for 1 ns, and independent replicates were generated by applying random velocities at its outset. Unrestrained simulations were then performed on equilibrated systems using OpenMM (68,69). The NPT ensemble was maintained, and the Monte Carlo barostat in OpenMM was used to isotropically regulate pressure with box scaling attempted every 25 integration steps.

Trajectory analysis was performed using standard utilities in the CHARMM (55) and R programs (62). Analysis techniques specific to this work are described below.

Ion sampling around GQs

Ion occupancy maps were generated as described elsewhere (45, 46, 47). The system volume was divided into discrete, 1-Å3 volume elements (voxels). Replicates from each unique system were combined and analyzed, and the location of each ion was assigned to the nearest voxel. We set an isosurface cutoff of ≥1%, meaning a voxel must be occupied by an ion for at least 1% of the frames to be visualized in the map. A ≥1% cutoff is permissive and enables a general description of ion sampling in flexible systems. We have found that this cutoff is sufficient to resolve discrete volumes of preferential ion binding.

Ion counting via volume Jacobian-normalized radial distribution functions

To calculate the excess ion accumulation around the c-kit1 GQ, we computed radial distribution functions, g(r), that were normalized by the volume Jacobian, J(r), at increasing distance from the GQ surface:where V is the total volume of the box, N is the total number of a given ion type, and N(r) is number of ions at distance r from the GQ surface. J(r) is the volume of a shell equidistant from the GQ surface (61) and does not monotonically increase (Fig. S2). Therefore, assuming a spherical or cylindrical Jacobian would lead to overestimation of ion condensation near the GQ surface (61). With this technique, GQ-ion g(r) can be computed by defining the GQ-ion distance as the minimal distance between DNA and a specific ion, generating GQ-ion distance histograms for each frame in the combined trajectories, and normalizing these histograms by J(r) (Fig. S2; 61). This method has been employed in previous studies with the Drude FF (44,46,60), leading to qualitative and quantitative descriptions of ion distributions around nucleic acids. A complete description of the theory and calculation can be found elsewhere (61). From the normalized g(r), the number of excess ions of a given type, N(r), can be computed fromwhere r is the bulk density (ions Å−3) of the ion of interest.

Local ion coordination

To compare K+, Na+, and Li+ ion coordination in the tetrad core, we calculated ion-tetrad relative distances. The relative distance of each tetrad and bound ion was centered on the tetrad core center of mass (relative distance = 0 Å). By combining replicates and comparing relative distances across all systems, we established a geometric criterion for “coplanar” and “antiprismatic” ion coordination. The relative positions of the tetrads were consistent across all systems (tetrad 1 = −3.4 ± 0.1 Å, tetrad 2 = 0 ± 0.1 Å, and tetrad 3 = 3.4 ± 0.1 Å), whereas relative ion positions were slightly more variable (ion 1 = −1.4 ± 0.8 Å and ion 2 = 2.1 ± 0.8 Å). Using these values, we described an ion to be “coplanar” if it was 0.8 Å or closer to the average position of the nearest tetrad and “antiprismatic” otherwise. These cutoffs were used to determine the impact of local ion coordination on the interaction energy of the tetrad core.

Interaction energy

Using the cutoffs described above, we performed interaction energy analysis between core cations and tetrad guanines as well as core cations and the surrounding water to determine the impact of local ion coordination. In total, three different interaction energy analyses were performed. First, we calculated interaction energies between core ions and the nearest tetrad guanine when K+, Na+, and Li+ were participating in coplanar or antiprismatic coordination. For this calculation, a tetrad-ion interaction was considered coplanar or antiprismatic if that tetrad had one ion (core ion 1 or core ion 2) interacting within the cutoff, resulting in “coplanar” and “antiprismatic” values for all three tetrads. Second, we calculated the interaction energy between core ions and all guanine bases in the tetrad core when K+, Na+, and Li+ were participating in coplanar or antiprismatic coordination. For the second calculation, in which we considered the entire tetrad core, both ions (core ion 1 and core ion 2) needed to fall within the cutoff to be considered coplanar or antiprismatic. For this reason, tetrad 1–3 interaction energies do not necessarily sum to the tetrad core interaction energies. The third interaction energy analysis performed was between the two core ions and water around the c-kit1 GQ (≤9 Å) as well as core ions and the nearest water molecules (any water atom ≤3.5 Å from an ion to account for the disordered nature of a partial hydration shell that exists around coordinated, coplanar ions).

Results and discussion

We have previously shown that the c-kit1 GQ can be stably modeled for 500 ns with the Drude FF (45). In solutions of ∼150 mM KCl, we showed that the tetrad core and coordinated ions were rigid and stable, that structurally important interactions in the long propeller loop (e.g., Gua17•Ade19) were preserved, and that bulk ions bind to previously described high occupancy sites (along the propeller loop, between Gua10 and Gua21, and above tetrad 1) (45). Here, we have simulated the c-kit1 GQ in ∼150 mM of KCl, NaCl, LiCl, and mixed salt solution for the first time with the Drude-2017 polarizable FF to determine if these observations differ as a function of monovalent cation identity and to thermodynamically and kinetically characterize the cationic environment around the c-kit1 GQ.

c-kit1 GQ stability and sampling

We evaluated the influence of cation type on c-kit1 GQ conformational sampling via root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), and dihedral sampling. In all systems, RMSD values were higher in the in the presence of LiCl and NaCl than in simulations with KCl, agreeing with the known GQ stabilization preference (K+ > Na+ > Li+). These observations are clearest when considering core nucleotide and base RMSD (Tables S1 and S2). Per-nucleotide RMSD analysis showed that these deviations were largest in core guanines near the 5′-end (guanines 2–4, 6–8, and 10), primarily when considering backbone heavy atoms (Figs. S3 and S4). Conversely, the RMSF values of all core guanines were low (Figs. S5 and S6), suggesting the backbone of core guanines near the 5′-end adopted slightly different conformations than the experimental structure but remained rigid. These deviations did not lead to core ion expulsion or notable disruption of Hoogsteen hydrogen bonds in bound and bulk K+ and Na+ simulations. However, guanine base fluctuations led to transient bifurcated hydrogen bonding in the tetrad core. These interactions did not persist between specific G:G basepairs, rather they occurred sporadically throughout the core. This phenomenon was most common in bound and bulk Li+ simulations and bulk-only simulations (Fig. S7). In each case, core hydrogen bond disruption increased the RMSD of the entire GQ structure, suggesting that local deviation or destabilization can influence the overall structure of the GQ. These findings agree with our previous work on the c-kit2 GQ (47). Because the c-kit1 and c-kit2 GQs have distinct topologies, the sensitivity of both structures to the monovalent cations in solution suggests that GQs with different folds and sequences are influenced by cation type.

Ion sampling around the c-kit1 GQ

We sought to describe the K+, Na+, and Li+ ion atmosphere around the c-kit1 GQ both qualitatively and quantitatively. To do so, we first generated ion occupancy maps as described in our previous studies (45, 46, 47). The ion occupancy maps revealed discrete regions of ion sampling that were similar in simulations with bound and bulk ions (Fig. 1) and bulk ions only (Fig. S8). Systems initialized with only bulk ions consistently led to ion binding to the tetrad core that resembled core ion coordination in the bound and bulk systems and was irreversible on the simulation timeframe. Further, there are similarities among the different ion types in sampling around the c-kit1 GQ that agree with previous simulation (45,70) and experimental data (71), specifically with respect to high occupancy above tetrad 1, around the long propeller loop (near O6 atoms of Gua17 and Gua18), and along the phosphodiester backbone between Gua10 and Gua21. There are interesting differences between bulk K+, Na+, and Li+ ion binding. First, bulk K+ aligned above tetrad 1 more frequently (76% of simulation frames) than Na+ (47%) and Li+ (31%), but these binding events affected the tetrad guanine base dipole moments differently. In KCl and NaCl, the nucleobase dipole moments increased in response to bulk ion alignment but were insensitive to alignment of bulk Li+ (Fig. S9). Similarly, we computed the dipole moments of each nucleobase in solutions of KCl, NaCl, and LiCl and found that the dipole moments of pyrimidine bases increased in systems of LiCl, whereas dipole moments of purine bases decreased (Table S3). Together, these results suggest that although Li+ systematically depolarizes guanine bases, Li+ binding above tetrad 1 does not perturb the electronic properties of guanine bases or purine bases in general, which could influence its ability to stabilize the native state of the c-kit1 GQ.

Figure 1

Ion interaction maps around the c-kit1 GQ, starting with bound and bulk ions in KCl (A) (gold), NaCl (B) (silver), and LiCl (C) (cyan). The isosurface value for ion sampling was set at an occupancy threshold of ≥1%. The displayed percentages reflect the persistence of each ion at the indicated location across the three replicate simulations, expressed as the fraction of snapshots in which an ion was aligned with the tetrad stem (see Methods). To see this figure in color, go online.

Li+ interaction maps (Fig. 1 C; Fig. S8 C) show more high occupancy sites than Na+ and K+, with sampling particularly enriched around backbone and sugar moieties (Fig. S10). These differences are interesting because external (noncore) ion-DNA interactions are known to contribute to GQ folding and stability of the folded state because cation binding to the negatively charged backbone should reduce electrostatic repulsion (8,9). However, structural deviations in our simulations were greatest in systems of LiCl. We have reported increased backbone interactions and structural deviation in simulations of the c-kit2 GQ (47) and show here that these observations extend to the c-kit1 GQ as well. It is possible that Li+ ions bind nonspecifically to the GQ surface and such nonspecific binding fails to stabilize important interactions in the folded GQ, but future simulations of the unfolded state would be required to test this hypothesis. Still, these results provide additional evidence that Li+ ions are not suitable for maintaining GQ stability because of nonspecific, diffuse backbone binding.

To better characterize the cation distributions around the c-kit1 GQ, we computed volume Jacobian-normalized g(r) in the extended mixed K+/Na+ and K+/Li+ systems. These calculations showed preferential accumulation of K+ ions over both Na+ and Li+, with Nex(K+) = 2.5 vs. Nex(Na+) = 2.0 in the mixed K+/Na+ system (Fig. 2 A) and Nex(K+) = 1.8 vs. Nex(Li+) = 1.7 in the mixed K+/Li+ system (Fig. 2 B). Thus, K+ can out-compete Na+ and Li+ for the surface of the c-kit1 GQ. The ion distributions of K+ and Na+ were similar across our simulations (Fig. 1, A and B; Fig. S8, A and B), suggesting that they compete for similar interactions. As noted above, the Li+ distributions were more diffuse (Fig. 1 C; Fig. S8 C), suggesting that the reduced difference in ion accumulation arises from the nonspecificity of Li+ binding to the c-kit1 GQ. The g(r) analysis also revealed that Na+ and Li+ ions were found closer to the GQ surface than K+ ions, and the difference was more pronounced in the case of Li+. We attribute this difference to the smaller radii of these ions, such that they make closer contact with GQ moieties.

Figure 2

Volume Jacobian-normalized radial distribution functions, g(r), of ions around the c-kit1 GQs in the (A) extended mixed K+/Na+ and (B) extended mixed K+/Li+ systems. The plotted g(r) for each solution condition is the average of three independent simulations of each system. To see this figure in color, go online.

Ion coordination in the tetrad core

To further understand the differences between K+, Na+, and Li+ ions in the c-kit1 GQ, we analyzed the local coordination of cations in the tetrad core via ion-tetrad relative distance distributions (Fig. 3). In KCl, the core ions maintained the expected bipyramidal antiprismatic coordination over 99% of the time. Core Na+ ions adopted both bipyramidal antiprismatic (40%) and coplanar coordination (60%), and core Li+ ions were coplanar with guanine bases over 99% of the time. This cation coordination agrees with outcomes of a previous density functional study (36), which indicates that quantum mechanical observations are upheld in a dynamic, solvated system and demonstrates the quality of the Drude FF for studying GQs. Although K+ ions remained in the canonical coordination geometry, Na+ ions moved down the tetrad axis to interact primarily with tetrads 2 and 3, whereas Li+ ions diffused either up (to interact with tetrads 1 and 2) or down (to interact with tetrads 2 and 3). As a result, the Li+ ion-tetrad relative distance distributions shown in Fig. 3 appear to have four separate traces, reflecting bimodal distributions for each of the core ions. Additionally, core Li+ ions were not always equidistant from the four coordinating O6 guanine atoms in the tetrad. Instead, they often interacted closer to two or three guanine bases in a tetrad such that the vector connecting the two ions was not coincident with the tetrad axis (Fig. 3). As a result of this asymmetry, the guanine bases furthest from a coplanar ion twisted to interact more closely with a coplanar Li+ in the neighboring tetrad, perturbing hydrogen bonding in the core (Fig. S7) and distorting nearby loop regions (Fig. S11). These results emphasize that coordination of ions in the tetrad core can influence the overall structure of the c-kit1 GQ and explains the elevated RMSD in the LiCl systems noted above.

Figure 3

Ion-tetrad relative distance distributions and core ion clusters show the position of the bound ions in the tetrad core. All distributions are centered on the tetrad core center of mass (relative distance = 0 Å). Ion-tetrad relative distances for all replicates were combined to produce distributions for KCl, NaCl, and LiCl. Cartoon renderings (right) illustrate dominant modes of ion coordination in these systems. To see this figure in color, go online.

Tetrad core interaction energy

Having found structural differences in the coordination of the K+, Na+, and Li+, we sought to determine if these changes affected the magnitude of nucleobase-ion interaction energy (EINT), which is expected to be the dominant stabilizing force in the native folded state of a GQ (23). Because Na+ and Li+ ions often adopted coplanar coordination with tetrad guanine bases and appeared to have a strong attraction to core guanine O6 atoms, we calculated core ion interaction energies when the ions were coplanar or bipyramidal antiprismatic in each tetrad and the entire core (Table 2). These values show that tetrad core ion interaction energies are stronger as a function of increasing charge density (EINT(Li+) < EINT(Na+) < EINT(K+)), which is a similar outcome as in our previous study on the c-kit2 GQ (47). Interestingly, tetrad-K+ interaction energies were strongest when the ion aligned in the canonical bipyramidal antiprismatic coordination, whereas for Na+ and Li+, tetrad-ion interaction energies were stronger when the ion was coplanar with coordinating tetrads. Although coplanar coordination strengthened individual tetrad-ion interaction energies, when considering the core as a whole, the tetrad core-ion interaction energies were strongest when ions were coordinated antiprismatically, regardless of ion type. Because tetrad core-ion interaction energies were weakened when ions aligned coplanar and Na+ and Li+ ions were frequently found with this alignment, these differences may contribute to the preference of c-kit1 for K+ ions in its native state.

Table 2

Interaction energies (EINT, kcal mol−1) between the core cations and the guanine bases of each tetrad in the c-kit1 GQ with KCl, NaCl, and LiCl

	Coplanar	Antiprismatic	ΔEINT
KCl
tetrad 1	−32.3 ± 4.1	−35.2 ± 3.8	2.9 ± 5.6
tetrad 2	−64.4 ± 3.8	−65.7 ± 4.2	1.3 ± 5.7
tetrad 3	−25.0 ± 5.2	−28.7 ± 4.2	3.7 ± 6.7
tetrad core	−112.8 ± 14.2	−129.5 ± 7.9	16.7 ± 16.2
NaCl
tetrad 1	−51.6 ± 8.0	−49.8 ± 8.9	−1.8 ± 12.0
tetrad 2	−83.3 ± 10.7	−78.0 ± 12.6	−5.3 ± 16.5
tetrad 3	−61.0 ± 5.8	−54.7 ± 7.7	−6.3 ± 9.6
tetrad core	−172.0 ± 16.2	−184.8 ± 11.4	12.8 ± 19.8
LiCl
tetrad 1	−81.6 ± 18.8	−59.4 ± 17.7	−22.2 ± 25.8
tetrad 2	−100.9 ± 13.2	−95.7 ± 12.7	−5.2 ± 18.3
tetrad 3	−101.8 ± 8.5	−86.1 ± 14.4	−15.7 ± 16.7
tetrad core	−212.0 ± 28.4	−240.1 ± 13.3	28.1 ± 31.4

Interaction energies were computed separately for frames in which ions were coordinated coplanar and bipyramidal antiprismatic in each tetrad (tetrads 1–3) as well as entire tetrad core (described in Methods). ΔEINT is calculated as EINT(coplanar) – EINT(antiprismatic).

To better understand why coplanar coordination of Na+ and Li+ ions to the tetrad persisted at the expense of weakened core interaction energies, we also calculated the interaction energies between bound ions and the water in the system (Table S4). This analysis revealed that ion-water interaction energies were similar across ion types when adopting antiprismatic coordination, which were weakly unfavorable, suggesting that the dominant stabilizing force in this geometry arises from GQ-ion interactions. Whereas K+-water interactions became more repulsive in the few frames corresponding to coplanar coordination, Na+- and Li+-water interaction energies became stronger, suggesting that favorable attraction from water partially enables this change in ion coordination. We then considered interactions between core ions and water molecules closest to the core ions (any atom within 3.5 Å of an ion). In this analysis, K+ ion-water interaction energies were repulsive and equal in magnitude in both coplanar and antiprismatic coordination (Table S5), suggesting that K+ ion-water interactions are not favorable in the GQ core regardless of local coordination. In contrast, Na+ ion-water interaction energy was stronger than K+ and stronger in antiprismatic coordination. This observation suggests that although interactions with the nearest water molecules favor antiprismatic coordination, longer-range water interactions (beyond 3.5 Å, or roughly the first hydration shell) can influence coordination preference and may help stabilize coplanar coordination of Na+. Last, Li+-water interaction energy was stronger than K+ and Na+ and stronger when in coplanar coordination, with the nearest water molecules contributing most strongly (Table S5). Together, these results suggest that attraction to surrounding waters can contribute to coplanar coordination in the tetrad core. Further, if Na+ and Li+ ions favor interactions with bulk water and coplanar coordination at the expense of weakened core interaction energies, it could also hinder ion binding to the tetrad core. These results challenge the “optimal fit” theory of K+ selectivity (29,30,36), as fit to the core also seems to be influenced by solvation properties.

Bulk cation binding and dehydration

To better understand the relationship of core ion alignment and hydration structure on K+, Na+, and Li+ ion binding, we analyzed simulations that were initialized without bound ions in the core. Ion binding occurred in all the simulations; however, we observed differences in the mechanisms by which K+, Na+, and Li+ bound to, and partitioned into, the tetrad core (Figure 4, Figure 5, Figure 6; Figs. S12–S17). In KCl systems, 1) the binding ions shed water molecules from their first solvation shell as they approached the open face above tetrad 1 and interacted with Thy12, 2) partitioned into the tetrad core and displaced water within the core, and 3) coordinated with the expected bipyramidal antiprismatic coordination (Fig. 4; Video S1). This stepwise process was observed in all replicates (Figs. S12 and S13; Videos S2 and S3). In systems of KCl, the average time for the first ion to bind was 4 ± 2 ns, resulting in two to three K+ ions coordinated to the tetrad core (Fig. S18).

Figure 4

K+ ion binding to the tetrad core. Trajectory snapshots were taken from a representative simulation in KCl. Water molecules colored red belong to the first hydration shell of the K+ ion that is in the act of binding, whereas pink water molecules initially hydrated the tetrad core. Time (ns) is denoted in the bottom right corner of each frame. To see this figure in color, go online.

Figure 5

Na+ ion binding to the tetrad core. Trajectory snapshots were taken from a representative simulation in NaCl. Water molecules colored red belong to the first hydration shell of the Na+ ion that is in the act of binding, whereas pink water molecules initially hydrated the tetrad core. Time (ns) is denoted in the bottom right corner of each frame. To see this figure in color, go online.

Figure 6

Li+ ion binding to the tetrad core. Trajectory snapshots were taken from a representative simulation in LiCl. Water molecules colored red belong to the first hydration shell of the Li+ ion that is in the act of binding, whereas pink water molecules initially hydrated the tetrad core. Time (ns) is denoted in the bottom right corner of each frame. To see this figure in color, go online.

Video S1. K+ binding and coordination in the tetrad core in Replicate 1

Relevant bulk ion binding and coordination events depicted between 0-2 ns, 25-27 ns, 121-123 ns, 158-160 ns, 164-166 ns, and 752-754 ns.

Video S2. K+ binding and coordination in the tetrad core in Replicate 2

Relevant bulk ion binding and coordination events depicted between 0-2 ns, 30-32 ns, 193-195 ns, 731-733 ns, 738-740 ns, and 791-793 ns.

Video S3. K+ binding and coordination in the tetrad core in Replicate 3

Relevant bulk ion binding and coordination events depicted between 8-10 ns and 36-38 ns.

The binding and partitioning of Na+ and Li+ ions differed from K+ ions. In NaCl systems, the binding ion did not always dehydrate before entering the core. Instead, it tended to retain 1–3 coordinated waters and pull these waters into the core while binding (Fig. 5; Fig. S14, S15, and S18; Videos S4, S5, and S6). Then, during the second or third Na+ binding event, the binding ion could shift perpendicular to the tetrad axis, allowing coordinated waters to exit the core (Video S4). Nevertheless, coordinated waters could be found between bound Na+ ions by the end of the 1-μs simulations. As described above and in our previous work with the c-kit1 GQ (45), Thy12 helped coordinate K+ binding in a clear, stepwise process. However, in NaCl solutions, binding ions could be coordinated by Ade1, Gua6, or Thy12. Further, in replicate 3, the tetrad core hydrogen bonding was disrupted early (Fig. S7), allowing ions to enter below tetrad 3 and through the long propeller loop. The average time for the first Na+ ion to bind was longer than K+ ions (24 ± 12 ns), resulting in 1–3 Na+ ions coordinated to the tetrad core (Fig. S18).

Video S4. Na+ binding and coordination in the tetrad core in Replicate 1

Relevant bulk ion binding and coordination events depicted between 8-10 ns and 42-44 ns.

Video S5. Na+ binding and coordination in the tetrad core in Replicate 2

Relevant bulk ion binding and coordination events depicted between 6-8 ns.

Video S6. Na+ binding and coordination in the tetrad core in Replicate 3

Relevant bulk ion binding and coordination events depicted between 40-42 ns and 56-58 ns.

These differences extended to systems of LiCl. Retention of first-shell waters upon Li+ binding to the c-kit1 GQ was even more exaggerated than in Na+ binding. Upon binding to the core, Li+ ions kept several coordinated waters and pulled them into the core (Videos S7, S8, and S9). In fact, bulk Li+ binding initially increased the number of water molecules in the tetrad core in all three replicates (Fig. 6 Fig. S15–S18). As in systems of NaCl, Li+ binding could be coordinated by Ade1, Gua6, or Thy12 and, in instances in which the tetrad core became disordered, Li+ ions entered below tetrad 3 and through the long propeller loop. Last, the average time for the first Li+ ion to bind was slower than both K+ and Na+ ions (54 ± 20 ns), resulting in two Li+ ions coordinated to the tetrad core (Fig. S18). The slow binding rate emphasizes the need for Li+ ions to overcome nonspecific interactions with other moieties in the c-kit1 GQ before they can partition into the tetrad core and potentially stabilize the structure.

Video S7. Li+ binding and coordination in the tetrad core in Replicate 1

Relevant bulk ion binding and coordination events depicted between 15-17 ns, 19-21 ns, 29-31 ns, and 37-39 ns.

Video S8. Li+ binding and coordination in the tetrad core in Replicate 2

Relevant bulk ion binding and coordination events depicted between 93-95 ns, 387-389 ns, and 466-468 ns.

Video S9. Li+ binding and coordination in the tetrad core in Replicate 3

Relevant bulk ion binding and coordination events depicted between 47-49 ns, 708-710 ns, and 711-713 ns.

In short, these single-salt simulations showed that 1) Na+ and Li+ can partition into the folded GQ without completely dehydrating, 2) bulk Na+ and Li+ binding mechanisms were more variable than that of K+, 3) more K+ binding events were observed than Na+ and Li+, and 4) bulk K+ binding was faster than either Na+ or Li+. It has long been believed that rate of ion dehydration influences GQ-ion preference and that cations have to be dehydrated completely to access the GQ core (34,72,73). Our simulations of the c-kit1 GQ showed that smaller Na+ and Li+ can retain water molecules and fit into the core; however, this phenomenon coincides with partial core disruption, and after binding, ions shift to let coordinated waters exit the core. Thus, dehydration of Na+ and Li+ appears critical for core integrity, but this dehydration can occur after the binding event. In addition, variable binding patterns in Na+ and Li+ may also be tied to rate of dehydration because these ions interacted near the tetrad core with their first solvation shell mostly intact for 1–13 ns before partitioning into the core, whereas these interactions with the open face occurred for <1 ns in all KCl replicates. These variable binding patterns might also suggest that Na+ and Li+ have a harder time forming specific contacts with the c-kit1 GQ. This idea is supported by analysis showing that K+ ions were more commonly found interacting with bases, Li+ ions more commonly interacted with backbone and sugar components, and Na+ ion interactions were somewhere between these trends (Fig. S10). Further, ion partitioning below tetrad 3 is counterintuitive because the long propeller loop partially blocks the tetrad core but was observed in previous simulations with the C36 FF (45). Because simulations with the C36 FF were more disordered than those with the Drude FF previously and simulations in NaCl and LiCl were more disordered than those with KCl here, it is possible that this ion binding trend is advantageous when the c-kit1 GQ becomes less stable, as might be the case when there are no core ions in systems of NaCl and LiCl. The facts that more ion binding events were observed in K+, followed by Na+ and Li+, and that first ion binding speeds were fastest for K+ than Na+, and Li+ are interesting and prompt further investigation of ion competition (discussed below). Importantly, each of these simulation outcomes support the stabilization trend K+ > Na+ > Li+.

Mixed salt binding in the c-kit1 GQ

Because GQs are unstable in the absence of ions, experimental analysis of ion binding and binding competition is difficult. To characterize the kinetics of cation binding and better understand if the differences described above regarding dehydration and ion binding speed contribute to preferential GQ:K+ binding, we performed several short simulations in solutions of mixed monovalent cations without bound core ions. In 40 mixed K+/Na+ simulations, 95 ion binding events were observed (64 K+ and 31 Na+ binding), and in 20 mixed K+/Li+ simulations, 44 ion binding events were observed (38 K+ and 6 Li+ binding). One-tail binomial tests (p < 0.05) performed using R software (62) showed that K+ ion binding was favored over both Na+ and Li+ (Fig. 7 A). The average speed of the first ion binding event was 15 ± 14 ns for K+ and 26 ± 24 ns for Na+ in mixed K+/Na+ simulations and 12 ± 10 ns for K+ and 34 ± 15 ns for Li+ in mixed K+/Li+ simulations (Fig. 7 B).

Figure 7

Bulk cation binding to the c-kit1 GQ. (A) Number of ion binding events in the mixed K+/Na+ system (left, n = 40, p = 4.62 × 10−4) and in the mixed K+/Li+ system (right, n = 20, p = 4.72 × 10−7). Asterisks denote the significance calculated from a one-tailed exact binomial test (α = 0.001). (B) Time of first ion binding event (ns) in the mixed K+/Na+ (K+, n = 23; Na+, n = 17) and K+/Li+ system (K+, n = 15; Li+, n = 5). x̄ denotes the average and SD of binding times. To see this figure in color, go online.

The observations that K+ preferentially binds to the tetrad core in solutions of K+/Na+ and K+/Li+ and that first ion binding event was fastest for K+ than Na+ and Li+ agree with our observations above and again support the stabilization trend of K+ > Na+ > Li+. Fewer and slower core ion binding events in Na+ and Li+ may be due, in part, to increased nonspecific interactions with the GQ surface. These outcomes also suggest a role for both kinetics, based on speed of dehydration, and thermodynamic preference at core binding sites, based on stronger core ion interaction energy and greater core integrity in antiprismatic ion coordination in preferential binding of K+ ions. Because we showed that Na+- and Li+-water interaction energies are strengthened in coplanar coordination and tetrad core disorder increases because of water in the core and coplanar coordination, preferred coordination of K+ over Na+ and Li+ in the c-kit1 GQ is likely dominated by the ability of each cation to dehydrate. Still, more work is required to understand if these observations extend to other GQ topologies, including the c-kit1 GQ with flanking nucleotides. We have shown previously that including non-GQ nucleotides can influence ion binding to the tetrad core and therefore could alter the ratio of binding events observed here (47). Similarly, further simulations are needed to fully understand how the relationship between hydration and coplanar coordination of ions influences the free energy profile of ion binding to the tetrad core.

Conclusions

Here, we have studied the effects of K+, Na+, and Li+ ions on the c-kit1 GQ by evaluating GQ structure, ion sampling, core energetics, ion dehydration and binding, and ion competition. From tracking structural deviation to comparing ion binding events in mixed salt solution, our simulations support the known GQ-ion specificity trend (K+ > Na+ > Li+) and provide new insights on preference for K+. As in our previous work with the c-kit2 GQ (47), we observed that K+ preferred antiprismatic coordination, Na+ ion shifted between antiprismatic and planar coordination, and Li+ adopted almost exclusively coplanar coordination. Interestingly, when Na+ and Li+ were coplanar, the tetrad-ion interaction energies and water-ion interaction energies were strengthened, whereas the entire tetrad core-ion interaction energies were weakened. This outcome suggests that rate of dehydration contributes to differences in local coordination of cations, which has previously been attributed to “optimal fit.” Further, our efforts in characterizing the ion atmosphere around the c-kit1 GQ showed that Na+ and Li+ interacted closer to the nucleic acid backbone and sugar than K+. Still, K+ ions preferentially accumulated around the GQ, suggesting that Na+ and Li+ may form nonspecific contacts with the GQ that influence the likelihood and speed of their binding as K+ ion binding was faster and more frequent than Na+ and Li+ in our simulations. In general, differences between K+, Na+, and Li+ simulations appeared to be strongly influenced by local ion coordination and specific versus nonspecific ion sampling, both of which correspond with ion dehydration properties. As such, we believe that preferred coordination of K+ over Na+ and Li+ in the c-kit1 GQ is dominated by the ability of the cations to dehydrate.