Comment on "Probing the thermodynamics of competitive ion binding using minimum energy structures".

We previously developed a theoretical framework to analyze ion selectivity of binding sites in macromolecules [?]. In a recent article [?], Rogers and Rempe stated in a note added proofs that this framework ignored an important free energy contribution [?]. A closer look shows that this statement is incorrect. To clarify the issue, let us consider a selective ion-binding site in a protein in equilibrium with bulk solvent. Assuming that ionic species i and j are present in solution, binding selectivity is governed by the relative free energy ΔΔGij=ΔGijsite−ΔGijbulk, where ΔGijbulk=[Gibulk−Gjbulk] is the free energy difference between ion i and j in the bulk solvent, and ΔGijsite=[Gisite−Gjsite] is the free energy difference between ion i and j in the binding site. The relative free energy of ion i and j in the binding site can be written as e−βΔGijsite=∫sitedXe−β[Ui(X)+ΔW(X)]∫sitedXe−β[Uj(X)+ΔW(X)] (1) where X represents the coordinates of a reduced subsystem comprising the ion and the ligands, Ui and Uj is the ionligand and ligand-ligand interactions for the system with ion types i and j, respectively, and ΔW(X) is an effective potential of mean force (PMF) that incorporates all the influence of the rest of the system (protein, membrane and solvent) on the local binding site subsystem. To an excellent approximation, ΔW does not depend on the ion type i and j if the charge of the ion does not change (i.e., no long-range effects). In the following, the framework is briefly summarized in a formulation that makes it easier to compare with that of ref [?]. First, we introduce the confining step-function Hc(X), which is equal to 1 for a subset of configuration and zero otherwise. The choice of the step-function will be specified later, but the following development is of general validity. We write the relative free energy of ions of type i and j in the binding site as, ΔGijsite=−kBTln[〈Hc〉j〈Hc〉i]+ΔGijr+ΔGijc (2) where the average of the ion-dependent confinement function is defined as, 〈Hc〉i=∫dXHce−β[Ui+ΔW]∫dXe−β[Ui+ΔW], (3) with a similar expression for 〈Hc〉j, and the relative free energies are, ΔGijc=−kBTln[∫dXHce−βUi∫dXHce−βUj], (4) and e−βΔGijr=[∫dXHce−β[Ui+ΔW]∫dXHce−βUi]×[∫dXHce−β[Uj+ΔW]∫dXHce−βUj]−1. (5) The significance of the various terms will be discussed below. The formulation based on Eq. (2) with various configurational restrictions is similar to previous treatments of ion solvation with constrained coordination numbers [?]. Such manipulations are common in computational methods for the calculations of absolute ligand affinities [?]. So far, the formulation is general, and a wide range of mathematical forms could possibly be used for the confinement step-function Hc(X). To pick a particular choice for Hc(X), we are motivated by the goal of clarifying the role of architectural forces on ion selectivity. For this purpose, we defined the architectural confinement for each atom k in the system as the smallest possible spherical volume Vk that encompasses all the dynamic excursions of that atom in the reference frame of the protein (no global translation/rotation) whether an ion of type i or j is bound in the site. In practice, the radius of the volumes Vk estimated from all-atom MD is typically on the order of 1.0-1.5 A [?]. This simple choice defines a minimal default model with probability distribution, that incorporates the generic effect of architectural confinement by the surrounding protein structure in an idealized fashion (without the protein, the volume Vk would be unbounded and the ligands would be allowed to flee away in space). By construction of the minimal model, 〈Hc〉j = 〈Hc〉i = 1, and ΔGijsite=ΔGijc+ΔGijr (6) The quantity ΔGijc given by Eq. (4) is the relative free energy governing ion selectivity arising in the subsystem confined by the function Hc(X), while ΔGijr accounts for the remaining (r) architectural forces from the surrounding protein. The latter may (or may not) serve to enforce a precise geometry to the atoms of the reduced subsystem. This specific choice Hc(X) allowed us to define the quantity ΔWc(X) ≡ –kBT ln[Hc(X)], which we have referred to as the “confinement” component of ΔW(X) in Eq. (1). The remainder, ΔWg ≡ ΔW – ΔWc, was referred to as the “geometry” component. This formulation provided a simple prescription to construct reasonable idealized generic reduced subsystems from the information extracted from all-atom MD simulations of complete biomolecular systems [?]. An important observation arising from this analysis is that selectivity can actually be realized via two extreme idealized mechanisms. The first case corresponds the snug-fit ideas familiar in host-guest chemistry in which the binding site provides a cavity of the suitable size to fit one ion but is unable to adapt to other ions. Selectivity is supported by architectural forces that enforce the local geometry of the site (ΔGijr is not negligible). The relative free energy ΔGij is roughly equal to the difference in the mean ionligand interaction energy. But the framework also helped explain how binding sites that are inherently flexible can also display robust ion selectivity [?]. We refer to this situation as the “confined micro-droplet” limit; ion selectivity emerges spontaneously at the level of the generic confinement without architectural forces imposing a precise geometry. In this case, ΔGijr is negligible and ΔGijc provides the dominant contribution to the relative free energy governing ion selectivity. The relative free energy of ion i and j is given, to a good approximation, by the difference in the mean ion-ligand and ligand-ligand interactions: ΔGijsite≈〈Ui〉i−〈Uj〉j. The trends in the confined micro-droplet limit can be robust and general; in an analysis of 1077 simplified reduced models comprising typical molecular groups, 39% displayed a non-trivial selectivity in the confined micro-droplet limit (i.e., |ΔΔGij| ≥ 2.0 kcal/mol) [?]. The theoretical formulation of Rogers and Rempe also relies on confinement step-functions, so-called indicator functions I(X; Cn) where Cn corresponds to the subset of allowed configurations of the system [?]. However, despite some similarities, there are also some important differences. In particular, the indicator functions introduced by the authors correspond to narrowly-defined ion-specific energy minimum coordination states. This choice is made to enable the effective use of the harmonic oscillator approximation of primitive quasichemical theory (pQCT) that is adopted by the authors to evaluate the I-constrained free energy ΔGijI (equivalent to Eq. (4)). In developing their formulation, the authors emphasize the importance of the free energy contribution arising from to the configuration constraint associated with the indicator functions, –kBT ln[〈I〉j/〈I〉i]. Indeed, the magnitude of the latter is likely to be significant because the averages 〈I〉i and 〈I〉j are expected to have different values and be much smaller than 1 in the case of narrowly-defined indicator functions. It is in this context that our previous analysis was suggested to be incorrect because it did not include a similar free energy contribution from confinement. In this regard, however, it seems that there was confusion. As shown above, the generic confinement Hc(X) in our analysis was specifically designed to encompass all the microscopic fluctuations of the binding site as observed in unbiased all-atom molecular dynamics simulations. This implies that the free energy contribution –kBT ln[〈Hc〉j/〈Hc〉i] is rigorously equal to zero by construction; see Eqs. (2)-(6). Thus, the assertion that the free energy contribution from confinement was ignored in our analysis is incorrect [?]. Ultimately, the goal of either of these approaches is to better evaluate the relative importance of ion-ligand and ligand-ligand interactions, Ui(X) and Uj(X), and architectural forces, ΔW(X), on the selectivity of ion binding sites in macromolecules. The framework formulated by Rogers and Rempe with narrowly-defined indicator functions o ers a valid approach to map reduced models onto well-defined energy minimum coordination states and enable the application of pQCT-like harmonic approximations [?]. Our approach was designed differently. We chose to first define a generic confinement from the fluctuations observed in all-atom MD simulations of membrane protein systems. This choice led to considerable simplifications, and helped clarify the construction of meaningful reduced subsystems from the information provided by the detailed MD simulations. The framework was used to analyze the underlying mechanism of ion selectivity in the S2 binding site of the KcsA K+ channel, the Na1 and Na2 binding sites of the LeuT leucine transporter [?]. More recently, it was used to clarify the effect of protonation on the K+ selectivity of the ion binding sites in the Na/K pump ATPase [?].