Strontium and barium in aqueous solution and a potassium channel binding site

Molecules that bind to potassium (K⁺) channels and block K⁺ ion permeation across cellular membranes can have detrimental or beneficial effects. Peptide toxins from several poisonous animals provide an example of the former.^1,2 For the latter case, drugs that block specific K⁺ channels hold promise for treating neurological disorders,³ autoimmune diseases,^3,4 and cancers.⁵ Alkaline earth metals also block K⁺ channels. As the simplest of the blocking molecules, divalent ions provide a favorable case for understanding the current issue of channel blocks. Two prominent ions used in K⁺ channel blocking studies are strontium (Sr²⁺) and barium (Ba²⁺). Comparable ion size and additional charge that lead to trapping compared to the permeant K⁺ ion have been proposed as an explanation for block by Ba²⁺.^6–8 This explanation naturally focuses attention on ion solvation by channel binding sites and free energies of ion transfer between water and those binding sites.

Despite close similarity of 0.16 Å in crystal radius and identical +2 charges,⁹ Ba²⁺ and Sr²⁺ ions exhibit different blocking behaviors in K⁺ channels. As the most potent of the metal ion blockers, Ba²⁺ binds preferentially to the innermost (S4) ion binding site in bacterial and eukaryotic K⁺ channels.^6,10 Ba²⁺ binding induces non-conducting events 10-1000-fold longer-lived than spontaneous channel closings, which allowed analysis of binding equilibria for K⁺ and Na⁺ ions.⁷ X-ray crystallographic studies of bacterial channels report that Ba²⁺ fully dehydrates upon transfer from aqueous solution and adopts the same octa-ligated conformation in the S4 binding site as permeant K⁺ ions (Fig. 1).^6,8,11,12 This similarity in solvation structure supports the proposed mechanism of divalent ion block.

By contrast, the slightly smaller Sr²⁺ ion binds less favorably to the S4 blocking sites of mammalian K⁺ channels^13,14 and does not block bacterial K⁺ channels.⁷ The difference in blocking behavior between Sr²⁺ and Ba²⁺ lacks an explanation and also undermines the traditional explanation for block. Could that difference be attributable to differences in hydration properties of these ions? Or could the difference instead be related to different solvation properties within the conserved S4 binding site of K⁺ channels?

Toward answering these questions, we previously reported the Ba²⁺ hydration structure and free energy.¹⁶ Here, we report a similar study of Sr²⁺ hydration. In addition, we carry out a direct comparison of Sr²⁺ and Ba²⁺ solvation properties at the innermost K⁺ channel blocking site, S4. To facilitate the comparison, we adapt the quasi-chemical free energy theory to describe solvation by the S4 binding site. In particular, we compute the absolute ion solvation free energies in water and the channel binding site to validate our approach and gain more insight over relative values computed more commonly. Our results predict neutral transfer free energies that suggest that +2 charge on ions similar in size to K⁺ does not lead to trapping behavior anticipated by earlier proposals. Our results also emphasize close similarities between the two divalent ions in terms of the solvation structure and free energy in the two distinct environments, as expected for ions of similar size and equivalent charge. This leads to the conclusion that differences in blocking behavior of Ba²⁺ and Sr²⁺ arise due to differences in dehydration and resolvation free energies, associated with the exchange of water ligands for channel ligands, instead of free energies of transfer from water to the binding site.

To probe the strontium hydration structure, we carried out ab initio molecular dynamics (AIMD) simulations and calculated radial distribution functions (RDF) between Sr²⁺ and water oxygens (O) (Fig. 2). The first shell is highly structured and well-defined, with eight (8) near-neighbor water molecules forming an inner hydration structure. A plateau on the running coordination number further confirms this tightly coordinated structure. The first minimum is spread over a broad range (3.4–3.8 Å), which is common for divalent¹⁷ and small monovalent ions.^18–20 

Compared to Ba²⁺, Sr²⁺ has a more narrow and well-defined first hydration shell.¹⁶ Also, the position of the first peak (2.63 Å) is slightly closer for Sr²⁺ than for hydrated Ba²⁺ (2.80 Å). Despite those small differences, eight (8) waters occupy the first hydration shells in both cases. This number differs from the first hydration shell of K⁺, which contains a total of 6 waters,^21,22 with only 4 tightly coordinated near neighbors that interact directly with the ion.²³ 

These simulation results compare well with the experimental results for Sr²⁺,^24–26 as was also the case for Ba²⁺ in earlier work.¹⁶ A slightly lower first shell occupancy of 7.5 obtained by a previous AIMD study²⁷ may be attributable to limited sampling from shorter simulation times (4 ps) compared with the current work (30 ps). The differences in RDF are consistent with the smaller ionic radius of Sr²⁺ compared to Ba²⁺, which pulls water molecules closer to the ion with higher charge density.

To evaluate hydration free energy, we applied quasi-chemical theory (QCT),^18,28,29 a complete statistical mechanical framework well-known for calculating such values.^30–32 This theory separates inner-shell interactions from outer-shells in a “divide and conquer” approach to obtain the solute excess chemical potential,

Detailed explanations of these terms and derivation of Eq. (1) are documented elsewhere.^28,33 In anticipation of ion solvation by ligands with multi-dentate modes of ion ligation (as in Fig. 1), we include a denticity factor m not given in earlier studies.

Briefly, the first term contains the equilibrium constant $Km*n(0)$ for association of n ligands with the ion to form an inner-shell cluster in an ideal gas, denoted by the superscript (0). An inner shell exists when ligands (e.g., waters) interact directly with the ion. If the ligand is multi-dentate, then m specifies the number of functional groups that coordinate with the ion (m = 1 for water). This association free energy term is a well-defined computational target that naturally takes into account full quantum chemical effects relevant to the chemical strength interactions between ions and ligands. Here we adopt the standard electronic structure software to evaluate this term and apply well-studied harmonic approximations^34,35 based on mechanically stable cluster structures. The harmonic approximation takes into account small atomic displacements from structures that represent minima on potential energy surfaces.

The density term, $ρH2Om*n$⁠, accounts for the availability of molecules (n waters with m = 1 functional groups, in this case) to act as ligands of the ion.

The second term gives the thermal probability, $pSr2+(m*n)$⁠, of observing an n-fold cluster in solution that ligates the ion with m functional groups. This probability can be estimated from AIMD simulation data. If only a single complex forms in solution, then the probability of observing that complex is 1 and the free energy contribution is zero. This situation can occur in strongly bound complexes formed from ligands interacting with multivalent ions such as Sr²⁺ or Ba²⁺. Alternatively, interactions between the cluster and surrounding environment may constrain the structure of the binding site. Structural constraints occur in many molecules, including ion channels.^36,37 While Eq. (1) can be written to account for any general structural constraint,¹⁸ here we only consider constraints on ligation number (m $*$ n). Since the value of the second term is negligible compared with the other terms, we neglect it here for computation of Sr²⁺ hydration free energy.

The third term accounts for the outer-shell solvation environment (e.g., here, liquid water) by taking the difference in excess chemical potentials between a cluster Sr$(H2O)m*n2+$ and n ligands. The clusters are treated as flexible molecular constituents of the solution. Evaluation of this term typically uses well-studied dielectric models, such as the polarizable continuum model (PCM) applied here.³⁸ This implicit solvation model is an important approximation, which is however balanced by the differences taken in the overall combination in Eq. (1).

The selection of an inner-shell occupancy and boundary are essential ingredients of QCT. For this hydration study, we explore a range of occupancies (n = 1–8) and boundary radii (λ = 2.4–3.6 Å) chosen from features in the radial distribution function, specifically the first hydration shell (Fig. 2). This elaborated study of the entire first hydration shell has not been done for Sr²⁺.

The application of Eq. (1) produces a range of ion hydration values (Fig. 3). The variation depends on the choice of occupancy and boundary used to define the inner hydration shell. Choosing a completely filled inner shell with n = 8 waters predicts a hydration free energy of $⟨μSr2+ex⟩$ = −331.8 kcal/mol that is independent of boundary radii and comparable to the experimental value of −333.6 kcal/mol.³⁹ The good agreement arises because a completely filled inner solvation complex reduces errors associated with the outer-shell term [Eq. (1)]. Specifically, a filled inner shell best satisfies the criteria for the successful application of an implicit solvent model of the outer-shell solvation environment.¹⁶ 

Decomposition of $μSr2+(ex)$ into inner and outer components (Fig. 3) illustrates that all terms of Eq. (1) naturally vary as a function of the inner-shell occupancy. Similar to observations of Li⁺,⁴⁰ Na⁺,⁴¹ K⁺,²³ Rb⁺,³⁰ Ba²⁺,¹⁶ and F⁻³² hydration by QCT analysis, the combined association free energy and ligand density term becomes more favorable as the number of ligands increases. Simultaneously, the balance of outer-shell contributions (difference between the cluster and n ligands) becomes less favorable with increasing n. The outer-shell contribution varies with the inner-shell boundary until the formation of the fully occupied inner hydration shell at n = 8 (Fig. 3). Then QCT, utilized with the approximations described above, produces an accurate ion hydration free energy that is also independent of the inner-shell boundary.

To summarize, Eq. (1) evaluated earlier with Ba$(H2O)1*82+$ and presently with Sr$(H2O)1*82+$ inner-shell clusters coupled to a polarizable continuum model of the surrounding environment produces accurate ion hydration free energies (Fig. 3).

Next, we consider ion hydration structures in the context of K⁺ channels. Earlier studies proposed that the octa-coordinated ion binding sites apparent in X-ray studies of K⁺ channels (Fig. 1) mimic the local hydration structure of permeant ions like K⁺ or Rb⁺.¹⁵ Studies of ion hydration by ab initio molecular dynamics showed that K⁺ and Rb⁺ local hydration structures are less than 8-fold coordinated.¹⁶ In other words, the K⁺ channel binding sites over-coordinate the permeant ions.^37,42 The binding sites mimic instead the hydration structure of Ba²⁺, a blocking ion. We find a similar result for the Sr²⁺ ion. Like Ba²⁺, Sr²⁺ forms stable 8-fold skewed cubic water clusters (Fig. 2).

From the point of view of drug design, the local ligand architecture of both divalent ions in water mimics binding sites for permeant ions in K⁺ channels (Fig. 1). The similarity in the local structure suggests that both divalent ions should fit similarly in potassium ion channel binding sites.

To examine ion solvation by K⁺ channels, we adapt the QCT analysis [Eq. (1)] for Ba²⁺ and Sr²⁺ solvation by the S4 binding site (Fig. 1). As a first task, we construct the explicit inner-shell solvation structures based on information from the K⁺ channel crystal structure with the S4 site occupied by Ba²⁺.¹¹ Specifically, the inner-shell solvation structure consists of four (n = 4) threonine (THR) amino acid residues that ligate the ion in the bidentate mode (m = 2). Using that information, we build inner-shell clusters Ba$(THR)2*42+$ and Sr$(THR)2*42+$⁠.

To validate these inner-shell structural models, we calculate the root-mean-squared deviation (RMSD) in distance of oxygen atoms between the THR residues of the crystal structure (1K4C) and the optimized Sr²⁺$(THR)2*4$ and Ba²⁺$(THR)2*4$ structures. All oxygen atoms move less than 1 Å after optimization of the electronic energy (Table I). This result confirms that the optimized structures provide a reasonable representation of the local structure for K⁺ ion channels occupied by Ba²⁺ in the S4 binding site (Fig. 4).

The inner-shell solvation structures represent simplified models of the ion channel. The study of simplified binding site models, with implicit or explicit inclusion of environmental effects, provides a natural route for learning about the underpinnings of selective ion binding under study here. Simplified models have been used in earlier studies of K⁺ vs. Na⁺ selectivity in similar models of interior K⁺ channel binding sites,³⁶ as well as binding sites formed of constrained water molecules.³⁵ In addition, the same model was used recently to study Ba²⁺ binding to variants of the S4 binding site.⁴³ Several new features differentiate the current work from those prior studies. First, we present for the first time the quasi-chemical theoretical framework for ligands with multi-dentate binding modes. Second, we account for the presence of liquid water in the outer-shell environment adjacent to the S4 binding site.

In the QCT analysis, the density term describes the availability of ligands to coordinate the ion. In analogy to the hydration studies, $ρTHR(2*4)$ may be taken as the density of pure liquid threonine at standard conditions. Liquid threonine density is the same order of magnitude as liquid water density (∼1 g/ml). Given that eight (8) oxygen ligands form a stable complex around Ba²⁺ both in water¹⁶ and in the S4 binding site,¹¹ the assignment of ligand density appears reasonable. A precise value for $ρTHR(2*4)$ is unimportant since the contribution to free energy depends on the natural log of density, $ln⁡ρTHR(2*4)$⁠.

Effects from the surrounding environment are taken into account in the probability and excess chemical potential terms. In the hydration studies, the calculated probability of observing the n = 8 complex is nearly 1 for both Sr²⁺ and Ba²⁺,¹⁶ making that term negligible in Eq. (1). Similarly, X-ray crystallography studies report a single structure, the n = 8 complex around Ba²⁺ in the S4 binding site,¹¹ suggesting that the probability term is also negligible in the channel environment for the 8-fold coordination case.

The excess chemical potential terms in the hydration problem modeled the aqueous environment as a polarizable continuum with a dielectric value of ϵ = 80. An aqueous environment borders the S4 binding site on the intracellular side in the form of a cavity with ∼20 waters in the closed state.¹⁵ In the open state, the cavity widens and becomes continuous with the bulk aqueous intracellular solution.⁴⁴ Toward the extracellular side, a single water borders the S4 binding site.¹⁵ 

A previous study of Ba²⁺ in the S4 binding site approximated the environment with ϵ = 1. Here instead we account approximately for the neighboring aqueous environment with ϵ = 80. The resulting ion solvation free energy is relatively insensitive to the precise value of ϵ for values between 80 and 30 (see the supplementary material). Furthermore, errors are balanced in taking the difference between excess chemical potentials of cluster and ligands. A more realistic model of the environment would treat the aqueous environment on one side of the S4 binding site differently from the other side, but this more complex model is reserved for future work.

The resulting QCT formula for ion (X) solvation by the S4 binding site describes 8 oxygen atoms from n = 4 THR residues that ligate the ions in the bi-dentate mode (m = 2),

Since the probability term p_X(2 $*$ 4) is likely negligible for the 8-fold coordination cases considered here, only the first term, the association free energy, and the last term, the outer-shell free energy, contribute to the ion solvation free energy in the S4 binding site, $μX−S4(ex)$⁠.

In the context of free energies for Sr²⁺ and Ba²⁺ transfer between water and K⁺ channel binding sites, the slightly smaller Sr²⁺ requires an additional ∼30 kcal/mol to take the ion out of the hydration environment (Fig. 5). The additional stability comes mainly from a stronger association free energy for Sr²⁺ in water relative to Ba²⁺, as expected due to the smaller size and higher charge density of Sr²⁺. This additional free energy must be compensated by the channel for both ions to have similar binding affinities.

Similar to the hydration results, Sr²⁺ is 30 kcal/mol more favorably solvated in the S4 binding site than Ba²⁺ (Fig. 5). Again, the additional stability comes from a stronger association free energy for the smaller Sr²⁺ in the S4 binding site. Similar outer-shell contributions are expected for the two ions in a given environment, as observed (Fig. 5), due to the similarity in cluster size, as well as identical ligand chemistry and number of ligands.

Comparing ion solvation in the S4 binding site relative to water (⁠$μX−S4(ex)$ versus $μX(ex)$⁠), the calculated transfer free energies show only minor differences between values of the same ion (Fig. 5). Interestingly, components of the free energy show substantial differences. The association free energy for both Sr²⁺ and Ba²⁺ in the S4 binding site of threonine oxygens is 30 kcal/mol more favorable than in water due to the stronger ion-ligand interactions from threonine oxygens. By contrast, the outer-shell contributions are ∼30 kcal/mol less favorable in S4 due to the larger size of the ion-ligand clusters in the S4 binding site. These opposing trends lead to a transfer free energy of ∼0 for both ions, suggesting that the blocking site environment can provide an additional 30 kcal/mol to stabilize the Sr²⁺ ion. Neutral transfer free energies also suggest that +2 charge on ions similar in size to the permeant K⁺ ion does not lead to trapping behavior (negative transfer free energies) anticipated by earlier hypotheses. A topic for future work consists of identifying the free energy barriers to divalent ion permeation within the channel, which likely arises due to under-coordination upon passage between the eight ligands of the S4 and S3 binding sites.

Strontium (Sr²⁺) and barium (Ba²⁺) ions are comparable in ion size and have identical +2 charge, but Ba²⁺ blocks bacterial potassium channels at the S4 binding site and Sr²⁺ does not. Our results suggest that in water, both ions adopt a local hydration structure consisting of 8-fold skewed cubic water clusters that mimic binding sites for permeant ions in potassium ion channels. But the smaller Sr²⁺ requires an additional ∼30 kcal/mol to take the ion out of the hydration environment relative to Ba²⁺. This additional free energy can be compensated by the S4 binding site, due to a more favorable association free energy term balanced by an equally less favorable outer-shell free energy, leading to a transfer free energy of ∼0. The transfer free energy for Ba²⁺ is also ∼0 and shows the same trends in contributions from the association and outer-shell terms. From the point of view of transfer free energies, these results suggest that both divalent ions should behave similarly and readily transfer between water and the S4 binding site of potassium channels. This view contrasts with earlier proposals that comparable size to the permeant K⁺ ion and +2 charge leads to blocking of K⁺ permeation by trapping of divalent ions.

The interesting observation on equivalent transfer free energies leads to the conclusion that the reason that Ba²⁺ blocks bacterial potassium channels in the S4 binding site, and Sr²⁺ does not, requires exploration of the free energies for the exchange of water ligands for channel ligands. The additional solvation free energy necessary to lose the same number of water molecules may hinder ligand exchange for Sr²⁺ compared to Ba²⁺. This conclusion is supported by the observation that both Sr²⁺ and Ba²⁺ block mammalian potassium channels, where the channel structure located in the intracellular environment may promote ligand exchange. This hypothesis will be tested in future work.

In summary, these results suggest that trapping of a divalent ion similar in size to the permeant K⁺ ion does not explain block of K⁺ permeation by Ba²⁺. Furthermore, differences in Sr²⁺ and Ba²⁺ blocking behavior do not arise due to the differences in free energies of transfer from water to the binding site. Instead, the free energies associated with the exchange of water ligands for channel ligands may account for the differences and will be the subject of future work.

The Vienna ab initio simulation package (VASP) was used with the PBE generalized gradient approximation⁵⁰ of electron density and the projector augmented wave (PAW) method⁴⁵ for core electrons. A constant volume (NVT) ensemble simulation was run for 30 ps. The first 20 ps of data were used for equilibration and the last 10 ps were used for analysis. The simulations were conducted at a temperature of T = 350 ± 21 K and ambient pressure, p = 1 atm, to be consistent with earlier studies of Ba²⁺ hydration.¹⁶ A cubic box of 12.417 Å dimension along each side, with 64 waters and a single Sr²⁺ ion, was simulated in periodic boundary conditions.

The system size and length of simulation time were sufficient to achieve a converged radial distribution function, confirmed by a separate classical MD simulation of the Ba²⁺ system simulated for 50 ns.¹⁶ In that box volume, the water density matches the experimental density of liquid water at the standard conditions of room temperature and pressure used in experimental structural studies. The boundaries contained a background charge to neutralize the overall charge of the system. With the PAW method, the valence electronic orbitals were expanded in plane waves with a high kinetic energy cutoff of 29.4 Ry (400 eV), and 10⁻⁴ eV was assigned as the convergence criteria for the electronic structure self-consistent iterations. A time step of 0.5 fs was used and a non-local van der Waals correlation function was applied to account for dispersion corrections using the vdw-DF2 method.⁴⁶ 

While experiments and simulation studies treat different ionic concentrations, experiments also show that the Sr²⁺ hydration structure is independent of ionic concentrations.⁴⁷ 

The initial configuration of the S4 blocking site, composed of four threonine residues $(Thr)2*4$⁠, was obtained from the crystal structure data for the bacterial potassium ion channel from Streptomyces lividans, KcsA (PDB ID: 1K4C).¹⁵ The K⁺ present in the crystal structure, K⁺$(THR)2*4$⁠, was replaced by either Ba²⁺ or Sr²⁺ ion. We selected the 1K4C structure because it has the highest resolution (2 Å) available for any potassium ion channel.

The KcsA crystal structure with Ba²⁺ at the S4 site is also available (PDB ID: 2ITD), but at lower resolution (2.5 Å).¹¹ Comparison between the 2ITD crystal structure and the computed lowest energy Ba²⁺$(THR)2*4$ structure indicates that the Ba²⁺–O distances in both structures are comparable (Fig. S2 of the supplementary material).

Electronic structure optimization and frequency calculations were done with the Gaussian suite of programs.⁴⁸ The TPSS functional and aug-cc-pvtz (O) and cc-pvdz (H) basis sets were used for hydration studies of Ba²⁺ and Sr²⁺. Corrections for electron correlation and electron delocalization were neglected here to facilitate comparison with earlier studies.^16,49 A smaller basis set (6-31+G*) was used for C, H, N, and O for binding site studies to decrease computational cost compared to the hydration calculations. To be consistent, we recalculated ion hydration free energy $μBa2+(ex)$ and $μSr2+(ex)$ using the 6-31+G* basis sets for H and O atoms of water (see Fig. 5 of the main text). Those values differ slightly compared with the higher accuracy results for Ba²⁺¹⁶ and Sr²⁺ hydration (Fig. S3 of the supplementary material) due to the smaller basis sets used. A polarizable continuum model (PCM) was used to represent solvent behavior with various dielectric constant (ε) values.

See supplementary material for comparison between the optimized structure of water and binding site clusters, μ^(ex) as a function of ε, comparison of 1ITD and optimized 1K4C clusters, and Gaussian input files.

We gratefully acknowledge support from the DTRA-JSTO CBD (IAA No. DTRA10027IA-3167) and Sandia’s LDRD program. Sandia National Laboratories (SNL) is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International, Inc. for the U.S. Department of Energy’s National Nuclear Security Administration under Contract No. DE-NA0003525. This work was performed, in part, at the Center for Integrated Nanotechnologies (CINT), an Office of Science User Facility operated for the U.S. DOE’s Office of Science by Los Alamos National Laboratory (Contract No. DE-AC52-06NA25296) and SNL. The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government.