Unveiling Zwitterionization of Glycine in the Microhydration Limit.

Charge separation under solvation stress conditions is a fundamental process that comes in many forms in doped water clusters. Yet, the mechanism of intramolecular charge separation, where constraints due to the molecular structure might be intricately tied to restricted solvation structures, remains largely unexplored. Microhydrated amino acids are such paradigmatic molecules. Ab initio simulations are carried out at 300 K in the frameworks of metadynamics sampling and thermodynamic integration to map the thermal mechanisms of zwitterionization using Gly(H2O) n with n = 4 and 10. In both cases, a similar water-mediated proton transfer chain mechanism is observed; yet, detailed analyses of thermodynamics and kinetics demonstrate that the charge-separated zwitterion is the preferred species only for n = 10 mainly due to kinetic stabilization. Structural analyses disclose that bifurcated H-bonded water bridges, connecting the cationic and anionic sites in the fluctuating microhydration network at room temperature, are enhanced in the transition-state ensemble exclusively for n = 10 and become overwhelmingly abundant in the stable zwitterion. The findings offer potential insights into charge separation under solvation stress conditions beyond the present example.

Introduction

Charge separation induced by water is among the fundamental processes in chemistry. A wealth of such processes is known, including but not limited to solvation of detached electrons and ion–electron pairs,[1−4] dissolution of salts,[5−7] self-ionization of water into the solvated proton and hydroxide,[8−12] or dissociation of acids,[13−15] which remain research topics of much current scientific interest. Common to all of them is the crucial role solvation plays foremost in the formation, but subsequently also in the stabilization of the charge-separated species by water. Both formation and stabilization of charges are advantageous if plenty of solvating water molecules are available, notably in bulk aqueous solution. The situation is radically different under the limit of having only a very few water molecules available: A piece of metallic sodium is known to most vigorously dissolve in a bucket of water by dissociating into Na+(aq) and e–(aq), while Na(H2O) are quite stable molecular complexes for sufficiently small cluster sizes n as shown by computation and experiment.[16−18] More recently, it has been shown that a minimum of four water molecules is required to dissociate a simple strong acid in the microhydration limit, HCl(H2O).[14,15,19−21]

While moving next to acids other than HCl is certainly appealing, a more complex situation is found in amino acids: They contain in one and the same molecule both the acidic proton-donating and the basic proton-receiving groups that are geometrically restrained by the molecular backbone. The well-known fact that amino acids exist in their neutral form in the gas phase, i.e., carrying charge-neutral −COOH and −NH2 groups, but readily convert into their zwitterionic form in bulk water supports the key role of solvation water in the charge-separation process leading to zwitterionization and thus to charged −COO– and −NH3+ groups. So far, extensive efforts have been made to understand zwitterion formation of both microhydrated and bulk solvated glycine as the role model for amino acids.[22−43] Whereas zwitterionization of glycine is well studied in the bulk phase limit[29] including its hydration properties therein[34,43] (e.g., the first hydration shell of glycine in water is shown to be constituted by about 7–8 water molecules around its charged groups[34,43]), not much is known in the microhydration limit on how zwitterionization can occur in terms of the underlying proton transfers and, second, which factors are contributing to zwitterion stabilization. To answer these questions requires not only the full mechanistic understanding of zwitterion formation in the microhydrated environment as such but also knowledge of the thermodynamics and kinetics by contrasting the generic scenario when the zwitterion is only metastable compared to its global stability. Despite long-standing efforts and most valuable computational insights into the “static geometry optimization limit ”,[22,24−28,33,35−38,41,42] these important questions still remain open even for the simplest possible amino acid at finite temperatures.

In what follows, we are going to address these questions based on free-energy surfaces (FESs) computed at room temperature to reveal the thermal mechanism and free energy of charge separation in glycine in the restricted solvation environment that is provided if only a few water molecules are available. To this end, we considered two microsolvated structures, Gly(H2O)4 and Gly(H2O)10, since it has been shown[38] that the undissociated (neutral) and charge-separated (zwitterionic) forms of glycine are their respective globally stable states, whereas the situation is much less clear-cut[38] for intermediate sizes 4 < n < 10 (as supported by new coupled cluster benchmark data; see SI Section 1 and Figure S1). Thus, contrasting these two limiting cases is expected to provide the most valuable insights into what finally favors zwitterionization in microhydrated Gly(H2O) clusters. To achieve this, we employed extensive ab initio molecular dynamics simulations[44] using the RPBE-D3 density functional (refer to SI Section 1 for validation of our methodology) together with metadynamics sampling and thermodynamic integration of the representative n = 4 and 10 clusters at 300 K (as detailed in the next section). We note at this early stage that the detailed picture of zwitterionization of microhydrated glycine thus obtained can be the basis for qualitatively rationalizing zwitterionization of other amino acids with neutral side chains, whereas such de/protonable functional groups will certainly add complexity. Moreover, our findings may not readily apply to biomacromolecules, such as polypeptides or proteins subject to solvation water restrictions, since the de/protonable terminal groups might be separated by large distances due to the presence of the backbone chain.

Results and Discussion

We first investigated the thermal reaction pathways by which neutral structures of Gly(H2O)4 and Gly(H2O)10 could zwitterionize at 300 K, given only 4 or 10 solvating water molecules. This is done by computing the free-energy landscapes (see the Methods and Models section for details) using the following two very general (dimensionless) collective variables (CVs): CV1 is defined as the difference of the coordination number (refer to the Methods and Models section for its definition) of nitrogen to all hydrogen atoms in the cluster, C[N–Hall], to the coordination number of O1, which lies in the trans position with respect to −NH2 and carries the hydrogen in the initial neutral cluster (see Figure for site labeling) again to all hydrogens present, C[O1–Hall], while CV2 is the coordination number of the other oxygen (O2) to all hydrogen atoms, C[O2–Hall]. We emphasize that this CV space is very general since it allows for proton transfers from/to any water molecule present in the cluster, thus not biasing transfers by selecting only a subset of water molecules, and at the same time involves all de/protonable sites of glycine. Furthermore, we note that the metadynamics simulations were exclusively performed to qualitatively explore the topology of the free-energy landscape of zwitterionization, whereas all reported quantitative free-energy differences (see below) were obtained from converged ab initio thermodynamic integration calculations (see the Methods and Models section) along suitable one-dimensional pathways as determined from metadynamics. Therefore, no further attempts were made to assess the convergence of these qualitative metadynamics simulations as the converged free-energy profiles, and thus the reported free-energy differences were finally obtained from thermodynamic integration calculations.

Figure 1

Free-energy surface at 300 K from ab initio metadynamics spanned by the generalized coordinates CV1 and CV2, see text, and representative configuration snapshots for neutral (N) ⇌ zwitterion (Z) interconversions for Gly(H20)4; see text. The blue arrows mark the locations of the shown structures on the FES; note that the neutral 4W-Noutcis and 4W-Nin conformers are indistinguishable in the (CV1, CV2) space, whereas the conformational difference between the 4W-Noutcis and 4W-Nouttrans isomers is visualized in Figure S3.

On performing the FES mapping first with 4W-Nout, a proton shuttling between the two carboxyl oxygen atoms via a water molecule was observed resulting in cis and trans isomers, 4W-Noutcis and 4W-Nouttrans, as indicated in Figure . These two isomers can be distinguished based on the location of protonated oxygen with respect to the −NH2 group, which is located in the cis position in 4W-Noutcis (corresponding to a dihedral torsional angle of θ(N–C–C–OH) ≈ 0°) and the trans position in 4W-Nouttrans (where θ ≈ 180°) as visualized using Newman projections in Figure S3. It is reassuring to sample this transformation since the existence of these particular cis and trans isomers of glycine has been found experimentally in the gas phase and also in the presence of one water molecule.[45−47] Our study, hence, shows that this isomerization can readily occur even in the presence of four microsolvating water molecules that clamp the carboxyl to the amino group. The neutral trans conformer, 4W-Nouttrans, is found to zwitterionize to yield 4W-Z via a water-mediated proton transfer from the carboxyl group to the amino group via three water molecules as depicted by the structure 4W-NZ. Interestingly, this metastable zwitterionic species 4W-Z can readily decay into another neutral species, 4W-Nin, via spontaneous direct proton transfer from its −NH3+ group to the neighboring oxygen of the carboxyl group, O2. This yields another conformer of neutral glycine in which the −COOH group points inward and donates an intramolecular H-bond to the amino group (see the 4W-Nin structure in Figure ); note that our metadynamics simulation does not resolve the neutral 4W-Nin conformer from 4W-Noutcis, which, however, does not impact on the free-energy profile of zwitterionization as reported below since that has been obtained using thermodynamic integration. The existence of this third neutral conformer was also experimentally observed[45,46,48] both for bare glycine in the gas phase as well as in Gly–water aggregates, which provides further confidence that our enhanced sampling approach faithfully maps the relevant conformational space.

In the case of 10W-Nout (see Figure ), on the contrary, no proton exchange between the two oxygens of the carboxyl group was observed, indicating the existence of only the trans conformer in larger clusters. In retrospect, this is expected as the nonprotonated O2 oxygen in 10W-Nout is involved in H-bonding with surrounding water molecules and, therefore, it has a reduced tendency to attract the proton from the neighboring oxygen atom. This is in contrast to 4W-Nout, where fewer water molecules around glycine limit the formation of such charge-stabilizing H-bonding interactions. The zwitterion, 10W-Z in Figure , is generated from 10W-Nout by a proton transfer cascade from the carboxyl oxygen to the amine nitrogen again via a chain of three water molecules as caught by the representative configuration 10W-NZ. Similar to what is observed in Gly(H2O)4, a proton transfer from −NH3+ to the carboxyl oxygen O2 interconverting 10W-Z into the neutral 10W-Nin species (see Figure ) is possible, yet only after surmounting a free-energy barrier.

Figure 2

Free-energy surface at 300 K from ab initio metadynamics spanned by the generalized coordinates CV1 and CV2, see text, and representative configuration snapshots for neutral (N) ⇌ zwitterion (Z) interconversions for Gly(H20)10; see text. The blue arrows mark the locations of the shown structures on the FES. The bifurcated H-bonded water bridge, see text, is depicted in 10W-Nouttrans, 10W-NZ, and 10W-Z by dashed green lines with the bifurcating H2O molecule enhanced, whereas the other H-bonds present in these structures are not shown for clarity. The H-bonds in 10W-Nin are shown using dashed red lines instead.

Given the similarities observed in the thermal reaction pathways for charge separation (i.e., Nout → Z) in Gly(H2O)4 and Gly(H2O)10 at 300 K, which are disclosed here for the first time, it is interesting to ask if the corresponding relative stabilities (i.e., thermodynamics) and energy barriers (i.e., kinetics) are also comparable. For this purpose, well-converged free-energy profiles have been computed for the Nout → Z zwitterionization of glycine using n = 4 and 10 water molecules by employing ab initio thermodynamic integration, the details of which can be found in SI Section 2. For n = 4, the thermal activation free energy for 4W-Nout → 4W-Z zwitterionization is about 12 kcal/mol, whereas the reverse barrier is only ≈4 kcal/mol (see Figure ). Thus, the neutral Gly(H2O)4 cluster is thermodynamically substantially more stable compared to its zwitterionic state. For n = 10, on the other hand, the forward and reverse barriers for the 10W-Nout → 10W-Z process are around 6 and 7 kcal/mol, respectively, implying that the zwitterionic state is thermodynamically stabilized by ≈1 kcal/mol (see Figure ); note that these thermal free-energy differences between neutral and zwitterionic species at 300 K are in qualitative agreement for both cluster sizes, n = 4 and 10, with the previously reported SCS-MP2 potential energy differences including zero-point vibrational corrections.[38] In conclusion, the charge-separated zwitterionic state is (marginally) stabilized for the larger microsolvated cluster, Gly(H2O)10, whereas the smaller cluster, Gly(H2O)4, prefers the neutral state at the level of thermodynamic stabilities at 300 K.

Figure 3

Free-energy profiles at 300 K from ab initio thermodynamic integration along the generalized coordinate CV1, see text and SI, for the 4W-Nout ⇌ 4W-Z (solid line with squares) and 10W-Nout ⇌ 10W-Z (dashed line with circles) neutral (N) ⇌ zwitterion (Z) interconversion processes. The left and right insets depict glycine itself in the corresponding neutral and zwitterionic states, respectively.

It now appears puzzling that the relative overall stabilities of the neutral and charge-separated states are different for Gly(H2O), given that the interconversion mechanisms are so similar for n = 4 and 10. Based on static optimizations of a vast ensemble of Gly(H2O) ground-state structures, it has been found[38] that a specific topological feature—namely, a bifurcated H-bonded water bridge, which, in addition to directly connecting the two functional groups, also connects the amino group to itself via a bifurcating water molecule—is required to render the zwitterion the global potential energy minimum (i.e., in the static 0 K limit and without providing any insights whatsoever into the charge-separation mechanism itself). At room temperature, however, H-bonded water networks are known to significantly fluctuate due to thermally driven H-bond breaking and making processes.[49] Thus, is there any hope that such bifurcated H-bonded water wires, being fragile features of H-bond networks, can be the answer? First of all, analyses of configurations underlying the free-energy surfaces show that the water molecules form an extended H-bond network that connects the two terminal groups in all cases (as visualized by the representative configurations in Figures and 2 and statistically analyzed in Figure S12). Second, however, only in the 10W-Nout to 10W-Z transformation process favoring the zwitterion is the presence of bifurcated H-bonded water bridges observed as shown in the snapshots in Figure . These bifurcated water wires are present in the n = 10 reactant state ensemble, slightly increase in number in the transition-state ensemble, and finally clearly dominate the H-bonding topology in the product state ensemble (see Figure S12 for quantitative analysis). The nontrivial finding is that this very topological feature of the microhydrating H-bond network is decisive for successful zwitterionization not only in the computational static zero Kelvin limit but also under realistic finite-temperature conditions where pronounced thermal fluctuations at 300 K induce enhanced H-bond dynamics.

For both cluster sizes, n = 4 and 10, the generation of nW-Nin from nW-Z opens up a new pathway for zwitterion to neutral transformation: Instead of nW-Z transforming into nW-Nout via the aforementioned water-mediated proton transfer chain, a direct proton transfer from −NH3+ to the O2 site of −COO– in nW-Z leads to another neutral conformation, nW-Nin. To determine the kinetic feasibility of this proton transfer pathway as well as the thermodynamic stability of nW-Nin with respect to nW-Z, we calculated the associated free-energy barriers for n = 4 and 10 as detailed in SI Section 3. The free-energy barrier for 4W-Z → 4W-Nin was obtained to be only ≈1 kcal/mol, while the reverse activation free energy is around 4 kcal/mol (see Figures and S7). Similarly, the forward and reverse barriers for 10W-Z ⇌ 10W-Nin interconversion are roughly 6 and 3 kcal/mol, respectively (see Figures and S7). However, 10W-Z is thermodynamically only marginally stabilized by ≈1 kcal/mol with reference to the neutral 10W-Nout species. Together, these results show that the zwitterionic state of Gly(H2O)10 is essentially only kinetically stabilized with respect to its neutral form. In contrast, the charge-separated 4W-Z state cannot be trapped at all since it is both thermodynamically and kinetically disfavored. These findings based on free-energy analyses were furthermore confirmed by running unbiased ab initio trajectory calculations of the 4W-Z and 10W-Z structures (refer to SI Section 3 for supporting data). Besides, we also scrutinized whether another interconversion channel is possible, involving the neutral Nin isomers instead (see SI Section 4). The computed free energetics, however, rule out kinetically favorable detour pathways for Z to Nout transformation via the Nin intermediate for both microhydration numbers (see Figures and S11). Overall, all of these analyses provide substantial evidence that the charge-separated zwitterionic form of Gly(H2O)10 is preferred under ambient gas-phase conditions mainly as a result of kinetic stabilization (since the lowest free-energy barrier for the neutral to zwitterion transformation is only 3 kcal/mol but twice as high for the lowest-energy reverse process), whereas the neutral form of Gly(H2O)4 is overwhelmingly populated at equilibrium since it is both thermodynamically and kinetically favored.

Figure 4

Free-energy barriers at 300 K from ab initio thermodynamic integration for conformational interconversions in (a) Gly(H2O)4 and (b) in Gly(H2O)10 in kcal/mol; note that kBT at 300 K corresponds to 0.6 kcal/mol. The thermodynamically stable states as well as the kinetically preferred pathways for charge recombination Z → N in (a) and for zwitterionization N → Z in (b) are highlighted in red.

Let us finally connect our findings on glycine microsolvation under ambient conditions to bulk solvation.[29] Unsurprisingly, water-mediated zwitterionization is also found in bulk solution, but now the charge-separated state is also thermodynamically stabilized by as much as about 10 kcal/mol (instead of being marginally stabilized by only ≈1 kcal/mol with n = 10 water molecules) together with an even slightly lower kinetic barrier of roughly 2 kcal/mol (versus ≈3 kcal/mol) to generate the zwitterion; recall that the thermal energy, kBT, is 0.6 kcal/mol at 300 K. This allows us to disclose the full sequence at room temperature from a situation where (i) zwitterionic glycine is both thermodynamically and kinetically unstable if too few water molecules are available (n = 4), to a scenario where (ii) the charge-separated state is mainly kinetically stabilized with more water molecules (n = 10), to the limit where (iii) the zwitterion is both thermodynamically and kinetically stabilized (n → ∞).

Conclusions and Outlook

In this work, we have determined the reaction pathways for intramolecular charge separation at finite temperature in two microsolvated amino acid clusters, namely, Gly(H2O)4 and Gly(H2O)10, through extensive ab initio simulations under ambient conditions. For both microhydration scenarios, the conversion of neutral into zwitterionic glycine takes place via a proton transfer cascade from the carboxyl to the amino group employing a short water wire. On comparing the free-energy barriers for all involved reactions as well as the relative free energies of the zwitterion with respect to neutral species, we demonstrate that the zwitterionic form of Gly(H2O)10 is thermodynamically only marginally stable at room temperature but kinetically stabilized, whereas it is the neutral form of the smaller cluster that is both thermodynamically and kinetically favored. The key feature allowing for sustained charge separation in the microhydration limit is revealed to be the formation of a specific topological feature of the H-bonding network. This is the establishment of bifurcated water bridges between the positively and negatively charged functional groups that are enhanced even at room temperature in the transition-state ensemble only for n = 10, and finally dominate over linear water wires in the fully charge-separated state. Bifurcated water bridges can only be built if a sufficient number of H-bonding water molecules is available beyond merely solvating the two charged groups by themselves. Overall, our study provides possible generic insights into zwitterionization of glycine and, as such, can be applied toward understanding charge-separation processes in other cases where both steric constraints and solvation stress are operational. Finally, it would be interesting to see in future studies how the presence of additional de/protonable groups in amino acids might impact the overall charge-separation mechanisms.

Methods and Models

System Setup and Computational Details

The minimum energy structures that correspond to the neutral and zwitterionic forms of Gly(H2O)4 and Gly(H2O)10, respectively, as reported in our previous work[38] were considered as the initial structures to launch the present ab initio molecular dynamics (AIMD) simulations.[44] The neutral and zwitterionic structures of Gly(H2O)4/Gly(H2O)10 will be denoted 4W-Nout/10W-Nout and 4W-Z/10W-Z, respectively. In addition, a minimum energy structure of one other conformer of the neutral form, where the −COOH group is rotated inward and forms a H-bond interaction with −NH2, was equilibrated and is denoted 4W-Nin/10W-Nin for Gly(H2O)4/Gly(H2O)10.

All simulations were performed using the CP2k suite of programs[50,51] that implements iterative Born–Oppenheimer propagation[44] to solve the electronic structure on the fly. The system was treated by density functional theory using the RPBE functional[52] supplemented by the D3 dispersion correction[53] (where the two-body terms together with zero damping are used) to account for the London (a.k.a. van der Waals) interactions. We chose the TZV2P triple-ζ Gaussian basis set with polarization functions[54] together with a plane wave cutoff of 500 Ry to represent the electron density, and the core electrons were described using norm-conserving separable dual-space Gaussian pseudopotentials.[55,56] A rather large simulation box of size (21 × 21 × 21) Å3 was considered in this study for all reported simulations, and the Martyna–Tuckerman Poisson solver[57] was used to apply finite-size rather than periodic boundary conditions. This is mandatory to fully suppress the long-range electrostatic interactions as required to properly compare free energies of zwitterionic versus the neutral Gly(H2O) clusters on equal footing. The temperature corresponding to 300 K was maintained using Nosé–Hoover chain thermostats.[58] A time step of 0.5 fs was chosen to integrate the equations of motion.

Free-Energy Sampling

The reaction pathways for the interconversion between different isomers were obtained by performing extended Lagrangian ab initio metadynamics simulations[59,60] using the same RPBE-D3 electronic structure method as that introduced in the previous subsection. Metadynamics is a powerful technique where sampling is enhanced in a space spanned by a set of predefined generalized coordinates, called the collective variables or CVs, which are assumed to provide a suitable coarse-grained description of the relevant free-energy subspace that describes all relevant rare events at the selected simulation temperature. To this end, a bias potential in terms of Gaussian hills is added to these CVs, which prevent the system from revisiting the already sampled points in the configurational space and by doing so allowing the system to visit unexplored regions including high (free)-energy regions corresponding to rare events. This eventually drives the system to escape the reactant minima via the lowest saddle point, enabling the reaction pathway to be determined. Additionally, an estimate of the underlying free-energy surface within the respective multidimensional CV subspace can be mapped out using the negative sum of the added biasing potentials. The interested reader can find more background and information on metadynamics sampling, for instance, in refs (44, 61−63).

We mainly used dimensionless CVs defined in terms of coordination numbers Cbetween a specific atom of species A and a set of other atoms J that belong to species B depending on their internuclear distance d[A – J]. Here, dAB0 is a fixed parameter that is defined based on the nature of interactions of the involved atoms, whereas p and q are integers that determine the steepness of the coordination number function. The value of dAB0 was chosen to be 1.3 Å in all simulations, whereas the values of both p and q were set to 12.

A Gaussian potential of height ≈1.0 kBT and width δ = 0.03 was used to fill the minima. To avoid the well-known “hill-surfing” problem in metadynamics,[64] we employed the adaptive time step approach, where a new Gaussian is added only after a displacement in CV space that exceeds 1.5 × δ with respect to the center of previously added Gaussian.[64] A mass of 50 amu was assigned to the CV, and a coupling parameter of 2.0 au was selected in the extended Lagrangian.[60] The CV temperature is maintained at 300 K via scaling with the accepted tolerance of ±200 K. These parameters have been validated in a set of previous studies where similar Grotthuss-like proton transfer reactions involving water molecules were scrutinized.[29,65−68]

Although metadynamics can be efficiently used for both accelerating rare events and reconstructing high-dimensional free-energy landscapes, the accuracy of the obtained free-energy difference and barrier estimates versus the underlying computational effort is a matter of concern.[61,62] An accurate description of barriers can only be achieved either by sampling many recrossing events involving all free-energy minima within the reaction subspace spanned by the respective CVs,[61,62] which is difficult or impossible to achieve on complex landscapes, or by performing well-tempered metadynamics where the growth rate of the bias potential is systematically decreased with simulation time,[69] where sampling becomes increasingly less efficient. In practice, both procedures are computationally demanding if the aim is to faithfully determine free-energy differences.

An alternative way to more easily obtain accurate estimates of free-energy differences along a one-dimensional transformation pathway only is to employ the traditional umbrella sampling[70] or thermodynamic integration (or potential of mean force)[71] methods along a predefined order parameter, which can be a typical CV, that connects the initial to the final state via a transition state. Therefore, in this work, instead of trying to quantitatively compute free-energy differences using metadynamics techniques, we employed ab initio thermodynamic integration[44] along a single CV. This was possible since analysis of the free-energy surfaces as obtained from ab initio metadynamics for both Gly(H2O)4 and Gly(H2O)10 revealed that orthogonal contributions to the interconversion of neutral and zwitterionic states as described by CV1 alone are very small according to Figures and 2 in the main text. For this purpose, multiple initial configurations along CV1 as given by the respective free-energy surfaces obtained from metadynamics were sampled; see SI Section 2 for definitions of CV1 and CV2. For each such configuration, a usual constrained ab initio simulation[44] at a fixed λ = CV1 value (of at least 50 ps length) was performed to provide a converged average ⟨f⟩λ of the force acting on the constraint for that particular constraint (after excluding the first 5 ps). The resulting one-dimensional free-energy profiles are calculated fromafter integrating the average constraint force from the respective initial to the final state.