Proton Transfer Pathways, Energy Landscape, and Kinetics in

Article pubs.acs.org/JPCB

Proton Transfer Pathways, Energy Landscape, and Kinetics in Creatine−Water Systems Olga Ivchenko,*,† Chris S. Whittleston,‡ Joanne M. Carr,‡ Petra Imhof,¶ Steffen Goerke,† Peter Bachert,† and David J. Wales*,‡ †

Department of Medical Physics in Radiology, Deutsches Krebsforschungszentrum (DKFZ, German Cancer Research Center), Im Neuenheimer Feld 280, 69120 Heidelberg, Germany ‡ University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, U.K. ¶ Institute of Theoretical Physics, Freie Universität Berlin, Arnimallee 14, 14195 Berlin, Germany ABSTRACT: We study the exchange processes of the metabolite creatine, which is present in both tumorous and normal tissues and has NH2 and NH groups that can transfer protons to water. Creatine produces chemical exchange saturation transfer (CEST) contrast in magnetic resonance imaging (MRI). The proton transfer pathway from zwitterionic creatine to water is examined using a kinetic transition network constructed from the discrete path sampling approach and an approximate quantum-chemical energy function, employing the self-consistent-charge density-functional tightbinding (SCC-DFTB) method. The resulting potential energy surface is visualized by constructing disconnectivity graphs. The energy landscape consists of two distinct regions corresponding to the zwitterionic creatine structures and deprotonated creatine. The activation energy that characterizes the proton transfer from the creatine NH2 group to water was determined from an Arrhenius fit of rate constants as a function of temperature, obtained from harmonic transition state theory. The result is in reasonable agreement with values obtained in water exchange spectroscopy (WEX) experiments.

■

(GAMT) deficiency.9,10 Furthermore, creatine analogues have proved to be anticancer agents that act synergistically with currently used chemotherapeutics.11 The purpose of this study is to understand the chemicalexchange process between creatine and water from computer simulation, and especially to investigate the mechanism of proton transfer and the associated kinetics. Creatine−water proton exchange can be characterized in terms of the activation energy required to transfer a proton from a creatine exchangeable group to a water molecule. This characterization is possible using geometry optimization techniques to identify transition states and pathways12,13 and, experimentally, by means of water exchange spectroscopy14 (WEX). The simulations of creatine in a water shell using an approximate treatment of the electronic structure yield insight into the basic mechanism of proton transfer occurring in metabolite−water systems and allow the corresponding kinetics to be quantified. We emphasize that our focus in the present work is on the structural transformations associated with proton transfer. We note that significant quantum dynamical effects are likely to play a role in determining the observed rate constants.15−21 Since the energetics of the rate-determining step are less well reproduced by the self-consistent-charge density-functional tight-binding (SCC-DFTB) potential, we do not consider such effects in the present work. Further characterization of the

INTRODUCTION Chemical exchange saturation transfer (CEST) enables molecular imaging in biological systems and thus provides a new magnetic resonance imaging (MRI) contrast for clinical diagnostics.42 The contrast is based on the chemical exchange between groups with “labile” protons in small metabolites (e.g., NH, NH2) and bulk water. Saturation-labeled proton transfer from the small solute pool to the much larger water pool produces an attenuation of the water signal intensity. We define ks→w as the corresponding chemical exchange rate for the transfer of a proton from the solute to water. The CEST effect depends on pH, temperature, and metabolite concentration. The advantage of the method lies in the imaging of specific metabolites (creatine,1 glucose,2 etc.) as well as the pH of tissue, which is typically altered in tumors.3,4 For human tumors, the average extracellular pH is pHe = 7.06 ± 0.05, while the average intracellular pH is pHi = 7.2.5,6 In comparison, a study using in vivo 1H-decoupled 31 P NMR spectroscopy in 17 healthy volunteers aged between 21 and 34 reported an average intracellular pHi of 7.11 ± 0.05 in normal calf muscle tissue,7 suggesting that tumors have a higher pHi.5 We consider metabolites of low molecular mass, which produce a CEST contrast in MRI, with creatine8 being the focus of the present work. Creatine possesses NH and NH2 groups, which can exchange protons with water. It plays a significant role in energy metabolism, and is seen at a depressed concentration in the tissues of patients with severe pathological conditions, such as guanidinoacetate methyltransferase © 2014 American Chemical Society

Received: March 21, 2013 Revised: January 18, 2014 Published: January 29, 2014 1969

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975

The Journal of Physical Chemistry B

Article

(Figure 2) using CHARMM.30,31 We chose the size of the shell in such a way that it contains enough water molecules to

most important kinetic pathways extracted from the SCCDFTB database may be possible using a more accurate treatment of the electronic structure, and we will report such calculations in the future. Accounting for tunnelling contributions to the rate constants will be more worthwhile within this framework. The NMR spectroscopic WEX method is the reverse of the CEST experiment and enables direct measurement of the exchange rate ks→w.22,23 The activation energy characterizing the proton transfer from the creatine guanidinium group to water was determined in the present work in simulations and compared to the experimental value of the effective activation energy from WEX.24 The term “effective” for the WEX activation energy is employed because this kinetic parameter is only valid for a certain buffer system (in this case, a phosphorus-based sodium−potassium buffer [phosphate buffered saline, PBS] was used in the WEX experiments). The present study probably represents the first attempt to model chemical exchange between creatine and water using a quantum mechanical treatment of the electronic structure.

■

Figure 2. End point minima used in the discrete path sampling procedure. Only the water molecules involved in the proton exchange are shown. The labile hydrogen is shown in yellow. Part A shows the zwitterionic structure, and part B shows deprotonated (negatively charged) creatine with an H5O2+ Zundel counterion.

Figure 1. Creatine in zwitterionic form encased by a water shell.

solvate the creatine but remained computationally tractable. Further solvation is not expected to be necessary, as the small barrier between solvation structures is washed out by the zeropoint motion of the proton,34 and further proton transfer, for example, between H3O+ and H5O2+, occurs on a femtosecond time scale and will not contribute to the observed proton transfer rate.35,36 Construction of a Kinetic Transition Network. Discrete paths37−41 (connected sequences of local minima and the intervening transition states) were obtained using the doubly nudged43 elastic band44−46 (DNEB) approach, followed by accurate refinement of transition states and a missingconnection algorithm,47 as implemented in OPTIM. DNEB is a chain-of-states method where a set of images is generated to interpolate between the initial and final minima of interest. The images are connected by springs, and “nudging” refers to a procedure for projecting out particular components of the true and spring gradients, which helps produce a better representation of the pathway.44−46 The doubly nudged implementation involves a further projection, which makes the band somewhat stiffer, and has been found to improve convergence of the highest energy images, which are the ones of primary interest.43 An initial linear interpolation was relaxed by gradient minimization for all the images simultaneously.43,48 Here a relatively large root-mean-square gradient threshold can be used as the stopping condition, since we only wish to identify the high-energy images that correspond to transition state candidates, and low-lying images in the vicinity of local minima are discarded.43,48 Any image corresponding to a local maximum was then considered as a starting point for accurate refinement to a transition state using the single-ended hybrid eigenvector-following approach.44,49,50 Here a transition state is defined geometrically as a stationary point with a single negative Hessian eigenvalue.51 This two-step procedure, involving DNEB and subsequent hybrid eigenvector-following refinement, has proved to be more efficient than converging all the images in a doubly ended search more tightly.43,48 The connections between pairs of minima defined by each new transition state were identified in terms of approximate

COMPUTATIONAL METHODS The simulations were performed using a geometry optimization framework, focusing on stationary points of the potential energy surface.12,13 Energies and first derivatives as a function of the nuclear coordinates were calculated from the selfconsistent-charge density-functional tight-binding (SCCDFTB) method25,26,28 using our interface between the OPTIM program29 and the CHARMM molecular simulation package.30,31 SCC-DFTB is an approximate quantum-chemical method derived from density functional theory (DFT) based on a second order expansion of the DFT total energy expression.32 This method was chosen, as it is 2−3 orders of magnitude faster than standard DFT calculations.33 Generation of Creatine Structure in Water. Creatine was surrounded by a water shell consisting of 35 water molecules (Figure 1), and initial end point structures were constructed in the zwitterionic8 and deprotonated forms

1970

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975

The Journal of Physical Chemistry B

Article

database,65−68 and through a pH range of 4−12.10,69 Our simulations therefore correspond to this regime.

steepest-descent paths, leaving the transition state parallel and antiparallel to the unique downhill direction. All the local minimizations, including the DNEB refinement, were performed with a modified version of the limited-memory Broyden−Fletcher−Goldfarb−Shanno (LBFGS) algorithm.52,53 Each new transition state and its corresponding minima were added to the database from which further pairs of minima were selected for connection attempts using the missing connection algorithm,47 where a distance-based metric is defined as input to the Dijkstra shortest path algorithm.54 The discrete path sampling (DPS) method37−41 for the simulation of rare events was used to construct a kinetic transition network (KTN)13,39,55−57 for the proton transfer process, as a connected set of minima and transition states. The resulting network was visualized using disconnectivity graphs,58,59 which provide an insightful representation of the barriers and structure of the energy landscape. The discrete path that makes the largest contribution to the DPS steadystate rate constant37 between the end points of interest was extracted using a graph-theoretical formulation60 and Dijkstra’s shortest-path algorithm.54 We refer to such discrete paths as the “fastest” paths, noting that this classification includes the conditional occupation probability of the reactant minimum.37 Global kinetic properties were extracted from the network using the NGT graph-transformation algorithm61 implemented in the PATHSAMPLE program.62 The experimental observation time scale will correspond to first order behavior if all other processes have relaxed. In terms of the master equation formulation of kinetics for a kinetic transition network, this scenario corresponds to a separation in time scales between the slowest process of interest and all other molecular scale relaxations.12,39 In the current work, this separation of time scales is clear-cut and the phenomenological rate constant, i.e., the experimental value, is understood as the rate of a single proton transfer from the NH2 group of zwitterionic creatine to any H2O in the water shell. This proton transfer rate is calculated between 0 and 37 °C and then used in an Arrhenius fit to determine the activation energy, Ea. In some cases, the NGT approach was applied in combination with a self-consistent regrouping procedure41 for the local minima. The regrouping converts the KTN in terms of stationary points of the potential energy surface into a network of free energy minima and transition states, where the members of each free energy minimum are assumed to achieve local equilibrium on a time scale that is much shorter than for transitions into an adjacent group, thus simplifying the kinetics in a physically justifiable manner. Choice of End Points. In the DPS approach, we have assumed that the proton transfer from zwitterionic creatine to water forming a local H5O2+ structure is the rate-determining step and that subsequent proton transfer to OH− is rapid in comparison, and can therefore be neglected. We believe this assumption to be reasonable, as the H5O2+ Zundel cation is the preferred structure at the air−water interface,34 with a 40% occurrence probability in bulk water,63 while the H3O+ hydronium ion is known to fluctuate between H5O2+ and H9O4+ complexes on a femtosecond time scale as a result of proton transfer.64 Further proton transfers from a particular water complex (e.g., H3O+, H5O2+, H7O3+, etc.) to OH− would be mediated by the water shell, producing a combinatorial number of routes for proton transfer from a specific water complex to OH−. The creatine zwitterionic structure is the dominant species at pH 7.3 according to the CHEBI

■

RESULTS AND DISCUSSION To characterize the mechanism of proton transfer for creatine in a water shell, we chose zwitterionic creatine and deprotonated creatine as two end points between which an initial path for proton transfer was calculated. We assigned these local minima as the reactant and product states; the structures are shown in Figure 2. The dynamics of proton transfer from the NH2 group of creatine to water occurs through a network of transition states. A connected path of minima and transition states between zwitterionic creatine and a deprotonated creatine structure was constructed.37−41 Once an initial connection between the two minima was obtained, the database of known structures was expanded using new connection attempts for minima of similar potential energy separated by high barriers. This procedure is designed to remove artificial frustration from the database caused by undersampling, as described in detail elsewhere.40 The resulting database, or kinetic transition network, contains 1839 minima and 2022 transition states, and was analyzed using the techniques described in the Computational Methods section. The phenomenological chemical exchange rate constant (obtained using the NGT algorithm61 and the self-consistent regrouping scheme41 with the regrouping energy threshold chosen to converge the value of the rate constant) is 7.5 × 10−6 s−1 at the physiological temperature T = 37 °C. When considering how to visualize the energy landscape, the distance between the nitrogen and hydrogen (dNH) of the creatine guanidinium group proved to be a useful order parameter for this proton-transfer reaction. The database was accordingly split into two parts, with zwitterionic structures defined by dNH < 1.5 Å and deprotonated creatine by dNH > 1.5 Å. The organization of the two sets of structures on the energy landscape was visualized by plotting a disconnectivity graph (Figure 3). Briefly, a disconnectivity graph has energy increasing up the vertical axis, and the position along the horizontal axis is arbitrary. The bottom of each vertical line corresponds to the energy of a particular minimum, and the minima are connected at the lowest energy for which they can interconvert, discretized at a series of threshold values to simplify the representation.58 The lines in the disconnectivity graph in Figure 3 have been colored according to whether they correlate with a zwitterionic structure (red) or deprotonated creatine (green). Figure 3 shows that the energy landscape of the creatine−water system clearly separates into two distinct regions according to the dNH order parameter, corresponding to zwitterionic (low energy; red in Figure 3) and deprotonated creatine (high energy; green in Figure 3). The clear separation of the disconnectivity graph into two sets of structures revealed in the figure shows that we have identified a good structural order parameter. It is important to note that this result is quite different from employing dNH as a “reaction coordinate”. The extraction of an effective first order rate constant, corresponding to the phenomenological rate constant observed experimentally, depends upon a separation of time scales for the relaxation of all phenomena relative to the process of interest. In the present case, the proton transfer will be slow compared to rearrangement time scales within the product and reactant regions identified by the disconnectivity graph using the structural order parameter. We note that dNH corresponds 1971

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975

The Journal of Physical Chemistry B

Article

creatine guanidinium proton is transferred from creatine to the free H5O2+ complex and thus to the water pool. With the two end points from this fastest individual path as reactant and product, the phenomenological rate constant at T = 37 °C is calculated as 7.5 × 10−6 s−1 with regrouping to convergence of the rate constant41 and 5.3 × 10−6 s−1 without regrouping. The former result is in quantitative agreement with the equivalent calculation using the original end points, as expected from the two-state nature of the energy landscape shown in the disconnectivity graph (Figure 3). Finally, taking all the zwitterionic minima from the database partition as reactants and all the remaining minima (deprotonated creatine) as products, the calculated rate constant at T = 37 °C without regrouping is 3.7 × 10−6 s−1. From the agreement to within an order of magnitude of all four calculated rate constants, we conclude that the chemical exchange rate constant for the SCCDFTB potential lies between 10−6 and 10−5 s−1 at the physiological temperature, and is not significantly affected by regrouping or the choice of end points. This value is much smaller than that found from experiment, as the WEX measurements were performed in PBS buffer rather than in pure water70 and barriers are known to be overestimated by the SCC-DFTB framework.32 Both of these factors are discussed below. The experimental chemical exchange rate at T = 37 °C obtained from CEST8 is ks→w = 950 ± 100 s−1, while WEX experiments24 yield ks→w = 1190 ± 119 s−1 at pH 7 and 1.2 ± 0.12 s−1 at pH 4. The activation energy Ea for proton transfer between zwitterionic and deprotonated creatine was determined by fitting rate constants in the temperature range 0−37 °C to the Arrhenius equation

k(T ) = A e−Ea / RT

(1)

with k(T) being the kinetic rate constant, T the absolute temperature, A a preexponential (or frequency) factor, and R the gas constant. The overall entropy of activation is included in the prefactor of the Arrhenius expression, while the activation energy in the exponential is an internal energy difference. Phenomenological rate constants k(T) were calculated between the original end points (see Figure 2), using the NGT algorithm61 and regrouping41 with the regrouping energy threshold chosen in each case to converge the value of the rate constant, for temperatures ranging from 0 to 37 °C. The temperature variation of k(T) in this range is accurately represented by an Arrhenius fit with an activation energy of EDFTB = 18.28 kcal/mol and preexponential factor of a ADFTB = 6 × 107 s−1. This activation energy corresponds to a proton transfer from the creatine guanidinium group to water in the model system consisting of creatine in a water shell of 35 water molecules. By assuming a dominant base-catalyzed and a weak buffercatalyzed proton exchange of the creatine guanidinium protons,22,24 WEX data allow us to determine the temperature dependence of ks→w, and hence to calculate the effective activation energy for proton transfer from zwitterionic creatine to a hydroxyl ion71 in the WEX experiment24 (see eq 14 in ref 24). The EDFTB value from simulations is about 2 times larger a than EWEX = 7.7 ± 1.77 kcal/mol obtained from the WEX a experiment.24 To assess this value of EWEX , we note that the a determination of the dependence of ks→w on pH and temperature in the WEX experiments is itself limited. There are two assumptions: (1) R1S (the longitudinal relaxation rate constant of the solute protons) does not depend on

Figure 3. Potential energy disconnectivity graph for creatine colored by protonation state. The graph includes 1839 minima and 2022 transition states that link them. Red minima correspond to creatine in the zwitterionic state and green minima to deprotonated creatine. An example structure from each protonation state is also shown. The disconnectivity analysis was performed in energy increments of δ = 1 kcal mol−1.

to an integrated path length, and is intrinsically multidimensional. The discrete path sampling approach avoids any need for projection onto lower dimensionality, and retains a faithful representation of all the barriers when calculating rates.37−39 The disconnectivity graph provides a direct visualization of how the structural order parameter successfully separates product and reactant states in this system, which is consistent with the observation of overall first order kinetics on the experimental observation time scale. Considering the database partitioned according to dNH as described above, the fastest individual path between any zwitterionic structure and any deprotonated creatine (without regrouping) makes the largest contribution to the rate constant for proton transfer in the creatine−water system. Figure 4 shows the corresponding potential energy as a function of integrated path length, and the pathway is shown in the inset. The transfer proceeds through the following steps: first, the proton from the creatine guanidinium group establishes a bond with the nearest water molecule, forming a hydronium ion, H3O+, and then, H3O+ combines with another water molecule to give an H5O2+ Zundel ion. Hence, the 1972

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975

The Journal of Physical Chemistry B

Article

Figure 4. The fastest proton transfer pathway identified for creatine as a function of relative potential energy, V. The zwitterionic starting point (A), transition state (TS), and charged final structure (B) are indicated.

temperature. In fact, the change of R1S with temperature is small compared to the change of ks→w with T. (2) The effect of buffer molecules in the solution on ks→w is negligible.72 In addition, the measurements were performed at low field (MR tomograph with B0 = 3 T), which required a long duration for the WATERGATE section73,74 of the WEX sequence, resulting in loss of signal. Given these assumptions, the accuracy of ks→w for creatine guanidinium proton transfer determined by WEX is about 10%.24 The experimental chemical exchange rates given above were measured for creatine-model solutions using phosphate buffered saline (PBS). A recent study of proton exchange in a serine protease70 has shown that using PBS increases proton exchange rates in comparison to the same system in water. As this effect is not currently accounted for in the WEX methodology, we expect that the proton transfer rate from the creatine NH2 group to water measured by WEX is overestimated with respect to the unbuffered solution. To assess the intrinsic error of the SCC-DFTB method, which results in the overestimated activation energy from the Arrhenius fit of the exchange rate (EDFTB = 18.28 kcal/mol) of a around 10 kcal/mol, another study has been undertaken for zwitterionic creatine solvated in a water shell consisting of 93 water molecules in a cubic box of side length 14 Å.27 The results will be discussed in detail elsewhere.27 We have determined the potential of mean force (PMF), which characterizes the height of the proton transfer barrier from the creatine guanidinium group to water using DFT umbrella sampling with the BLYP functional.27 The free energy profile of the proton transfer from the creatine guanidinium group to water was obtained from this PMF. The resulting proton transfer barrier of EaDFT = 3.012 ± 0.007 kcal/mol27 (12.60 ± 0.001 kJ/mol) is relatively small, and suggests that spontaneous proton transfer from the creatine guanidinium group to water will occur. This transfer was subsequently detected in an unbiased molecular dynamics (MD) run of 348 ps. From the ratio of the number of frames corresponding to the transferred state (deprotonated creatine) and the zwitterionic creatine from an unbiased MD run, a free energy difference of 11.3 kT (about 7 kcal/mol) can be estimated.27 The free energy difference between the two states, calculated

from the unbiased production MD with DFT and BLYP when the spontaneous proton transfer occurs, provides a separate lower bound to the barrier, which differs by 11 kcal/mol from the activation energy obtained from SCC-DFTB. Taking into account the typical error from the PMF in the umbrella sampling simulations with DFT and the BLYP functional, a variation of about 13 kcal/mol is observed between the two methods. The energy barrier obtained from the simulations, EDFTB , is a too large, and this discrepancy probably results from the SCCDFTB potential. The barrier calculated within the same framework for proton transfer from an sp3 nitrogen atom to either oxygen or nitrogen atoms is 10 kcal/mol above the correct value.32 The activation energy calculated for the cellobiohydrolase I glycosylation reaction using SCC-DFTB is 17.5 kcal/mol, which is also consistent with the present results.75 Nevertheless, structurally, SCC-DFTB has been found to describe hydrogen-bonded geometries to within 0.1 Å.32 It is also reported to reproduce the protonation and deprotonation energies of many compounds reasonably well, compared to values obtained using DFT with the B3LYP functional.32,76,77 For the present purposes, it is the proton transfer mechanism that is the focus of interest, and this pathway is likely to be better reproduced than the barrier height.

■

CONCLUSION The mechanism for proton transfer between the creatine guanidinium group and water has been investigated using computer simulations and an approximate quantum-chemical energy function, SCC-DFTB. A kinetic transition network was constructed using geometry optimization techniques, and the resulting potential energy surface was visualized using disconnectivity graphs. The graphs show a clear energetic separation between zwitterionic and deprotonated creatine structures. A detailed mechanism for proton transfer was obtained by analyzing the corresponding pathways. To characterize the global kinetics of creatine−water proton transfer in our simulations, we calculated the rate constant as a function of temperature, and obtained the activation energy from an Arrhenius fit. This activation energy, EaDFTB = 1973

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975

The Journal of Physical Chemistry B

Article

Polymer-Sensor: a Diagnostic Tool Applicable to Creatine Deficiency Syndrome. Biomed. Chromatogr. 2007, 21, 976−986. (11) Wyss, M.; Kaddurah-Daouk, R. Creatine and Creatinine Metabolism. Physiol. Rev. 2000, 80, 1108−1182. (12) Wales, D. J. Energy Landscapes; Cambridge University Press: Cambridge, U.K., 2003. (13) Wales, D. J. Energy Landscapes: Some New Horizons. Curr. Opin. Struct. Biol. 2010, 20, 3−10. (14) Mori, S.; Johnson, M. O.; Berg, J. M.; van Zijl, P. C. M. Water Exchange Filter (WEX Filter) for Nuclear Magnetic Resonance Studies of Macromolecules. J. Am. Chem. Soc. 1994, 116, 11982− 11984. (15) Cui, Q.; Karplus, M. Quantum Mechanics/Molecular Mechanics Studies of Triosephosphate Isomerase-Catalyzed Reactions: An Effect of Geometry and Tunneling on Proton-Transfer Rate Constants. J. Am. Chem. Soc. 2002, 124, 3093−3124. (16) Gao, J.; Truhlar, D. G. Quantum Mechanical Methods for Enzyme Kinetics. Annu. Rev. Phys. Chem. 2002, 53, 467−505. (17) Pu, J.; Gao, J.; Truhlar, D. G. Multidimensional Tunneling, Recrossing, and the Transmission Coefficient for Enzymatic Reactions. Chem. Rev. 2006, 106, 3140−3169. (18) Fleming, G. R.; Huelga, S. F.; Plenio, M. B. Focus on Quantum Effects and Noise in Biomolecules. New J. Phys. 2011, 13, 1−5. (19) J. P Bothma, J. B. G.; McKenzie, R. H. The Role of Quantum Effects in Proton Transfer Reactions in Enzymes: Quantum Tunneling in a Noisy Environment? New J. Phys. 2010, 12, 1−27. (20) Borgis, D.; Hynes, J. T. Dynamical Theory of Proton Tunneling Transfer Rates in Solution: General Formulation. Chem. Phys. 1993, 170, 315−346. (21) Truhlar, D. G. Tunneling in Enzymatic and Nonenzymatic Hydrogen Transfer Reactions. J. Phys. Org. Chem. 2010, 23, 660−676. (22) Zhou, J.; Wilson, D. A.; Sun, P. Z.; Klaus, J. A.; Van Zijl, P. C. M. Quantitative Description of Proton Exchange Processes Between Water and Endogenous and Exogenous Agents for WEX, CEST, and APT Experiments. Magn. Reson. Med. 2004, 51, 945−952. (23) Mori, S.; Abeygunawardana, C.; Van Zijl, P. C. M.; Berg, J. M. Water Exchange Filter with Improved Sensitivity (WEX II) to Study Solvent-Exchangeable Protons. Application to the Consensus Zinc Finger Peptide CP-1. J. Magn. Reson. 1996, 110, 96−101. (24) Goerke, S.; Zaiss; M.;.Bachert, P. Characterization of Creatine Guanidinium Proton Exchange by Water-Exchange (WEX) Spectroscopy for Absolute-pH CEST Imaging in vitro. NMR Biomed., DOI: 10.1002/nbm.3086. (25) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260−7268. (26) Frauenheim, T.; Seifert, G.; Elsterner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. A Self-Consistent Charge Density-Functional Based Tight-Binding Method for Predictive Materials Simulations in Physics, Chemistry and Biology. Phys. Status Solidi B 2000, 217, 41−62. (27) Ivchenko, O.; Bachert; P.; Imhof, P. Umbrella sampling of proton transfer in a creatine-water system. Manuscript in preparation. (28) Frauenheim, T.; Seifert, G.; Elstner, M.; Niehaus, T.; Köhler, C.; Amkreutz, M.; Sternberg, M.; Hajnal, Z.; Carlo, A. D.; Suhai, S. Atomistic Simulations of Complex Materials: Ground-State and Excited-State Properties. J. Phys.: Condens. Matter 2002, 14, 3015− 3047. (29) Wales, D. J. OPTIM: A Program for Optimising Geometries and Calculating Pathways. http://www-wales.ch.cam.ac.uk/OPTIM/. (30) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187−217. (31) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614.

18.28 kcal/mol, is higher than the result from experimental WEX data, EWEX = 7.7 kcal/mol, as anticipated from systematic a errors previously observed for the SCC-DFTB potential32 and the buffer effects70 not accounted for in the WEX methodology. Nevertheless, our results provide fundamental insights into the mechanism of proton transfer in metabolite−water systems and suggest a possible approach for the quantification of exchange rates in CEST.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was carried out within the HPC-EUROPA2 project (project number: 228398) with the support of the European Commission - Capacities Area - Research Infrastructure. It made use of the facilities of HECToR, the U.K.’s national highperformance computing service, which is provided by UoE HPCx Ltd at the University of Edinburgh, Cray Inc, and NAG Ltd, and funded by the Office of Science and Technology through EPSRC’s High End Computing Programme. We gratefully acknowledge the stipend for Olga Ivchenko from the International Max Planck Research School of Quantum Dynamics. This research was also funded by ERC grant RG59508 and EPSRC grant RG58073.

■

REFERENCES

(1) Kogan, F.; Haris, M.; Singh, A.; Cai, K.; Debrosse, C.; Nanga, R. P. R.; Hariharan, H.; Reddy, R. Method for High-resolution Imaging of Creatine in vivo Using Chemical Exchange Saturation Transfer. Magn. Reson. Med. 2014, 71, 164−172. (2) Chan, K. W. Y.; McMahon, M. T.; Kato, Y.; Liu, G.; Bulte, J. W. M.; Bhujwalla, Z. M.; Artemov, D.; van Zijl, P. C. M. Natural Dglucose as a Biodegradable MRI Contrast Agent for Detecting Cancer. Magn. Reson. Med. 2012, 68, 1764−1773. (3) Zhou, J.; Payen, J.-F.; Wilson, D. A.; Traystman, R. J.; Van Zijl, P. C. M. Using the Amide Proton Signals of Intracellular Proteins and Peptides To Detect pH Effects in MRI. Nat. Med. 2003, 9, 1085− 1090. (4) Garcia-Martin, M. L.; Martinez, G. V.; Raghunand, N.; Sherry, A. D.; Zhang, S.; Gillies, R. J. High Resolution pH(e) Imaging of Rat Glioma using pH-dependent Relaxivity. Magn. Reson. Med. 2006, 55, 309−315. (5) Swietach, P.; Vaughan-Jones, R. D.; Harris, A. L. Regulation of Tumor pH and the Role of Carbonic Anhydrase 9. Cancer Metastasis Rev. 2007, 26, 299−310. (6) Engin, K.; Leeper, D. B.; Cater, J. R.; Thistlethwaite, A. J.; Tupchong, L.; McFarlane, J. D. Extracellular pH Distribution in Human Tumours. Int. J. Hyperthermia 1995, 11, 211−216. (7) Bachert, P. Habilitation thesis, University of Heidelberg, 1995. (8) Haris, M.; Nanga, R. P. R.; Singh, A.; Cai, K.; Kogan, F.; Hariharan, H.; Reddy, R. Exchange Rates of Creatine Kinase Metabolites: Feasibility of Imaging Creatine by Chemical Exchange Saturation Transfer MRI. NMR Biomed. 2012, 25, 1305−1309. (9) Schulze, A.; Bachert, P.; Schlemmer, H.; Harting, I.; Polster, T.; Salomons, G.; Verhoeven, N.; Jakobs, C.; Fowler, B.; Hoffmann, G.; Mayatepek, E. Lack of Creatine in Muscle and Brain in an Adult with GAMT Deficiency. Ann. Neurol. 2003, 53, 248−251. (10) Sharma, P. S.; Lakshmi, D.; Prasad, B. B. Molecularly Imprinted Solid-Phase Extraction Combined with Molecularly Imprinted 1974

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975

The Journal of Physical Chemistry B

Article

(32) Elstner, M. The SCC-DFTB Method and Its Application to Biological Systems. Theor. Chem. Acc. 2006, 116, 316−325. (33) Hou, G.; Zhu, X.; Cui, Q. An Implicit Solvent Model for SCCDFTB with Charge-Dependent Radii. J. Chem. Theory Comput. 2010, 6, 2303−2314. (34) Vacha, R.; Buch, V.; Milet, A.; Devlin, J. P.; Jungwirth, P. Autoionization at the Surface of Neat Water: Is the Top Layer pH Neutral, Basic, or Acidic? Phys. Chem. Chem. Phys. 2007, 9, 4736− 4747. (35) Woutersen, S.; Bakker, H. J. Ultrafast Vibrational and Structural Dynamics of the Proton in Liquid Water. Phys. Rev. Lett. 2006, 96, 138305-1−138305-4. (36) Kirchner, B. Eigen or Zundel Ion: News from Calculated and Experimental Photoelectron Spectroscopy. ChemPhysChem 2007, 8, 41−43. (37) Wales, D. J. Discrete Path Sampling. Mol. Phys. 2002, 100, 3285−3306. (38) Wales, D. J. Some Further Applications of Discrete Path Sampling to Cluster Isomerization. Mol. Phys. 2004, 102, 891−908. (39) Wales, D. J. Energy Landscapes: Calculating Pathways and Rates. Int. Rev. Phys. Chem. 2006, 25, 237−282. (40) Strodel, B.; Whittleston, C. S.; Wales, D. J. Thermodynamics and Kinetics of Aggregation for the GNNQQNY Peptide. J. Am. Chem. Soc. 2007, 129, 16005−16014. (41) Carr, J. M.; Wales, D. J. Folding Pathways and Rates for the Three-Stranded beta-Sheet Peptide Beta3s using Discrete Path Sampling. J. Phys. Chem. B 2008, 112, 8760−8769. (42) Van Zijl, P. C. M.; Yadav, N. N. Chemical Exchange Saturation Transfer (CEST): What is in a Name and What Isn’t? Magn. Reson. Med. 2011, 65, 927−948. (43) Trygubenko, S. A.; Wales, D. J. A Doubly Nudged Elastic Band Method for Finding Transition States. J. Chem. Phys. 2004, 120, 2082− 2094. (44) Henkelman, G.; Jónsson, H. A Dimer Method for Finding Saddle Points on High Dimensional Potential Surfaces Using Only First Derivatives. J. Chem. Phys. 1999, 111, 7010−7022. (45) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method For Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (46) Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978−9985. (47) Carr, J. M.; Trygubenko, S. A.; Wales, D. J. Finding Pathways between Distant Local Minima. J. Chem. Phys. 2005, 122, 234903-1− 234903-7. (48) Sheppard, D.; Terrell, R.; Henkelman, G. Optimization Methods for Finding Minimum Energy Paths. J. Chem. Phys. 2008, 128, 134106-1−134106-10. (49) Munro, L. J.; Wales, D. J. Defect Migration in Crystalline Silicon. Phys. Rev. B 1999, 59, 3969−3980. (50) Kumeda, Y.; Munro, L. J.; Wales, D. J. Transition States and Rearrangement Mechanisms from Hybrid Eigenvector-Following and Density Functional Theory. Application to C10H10 and Defect Migration in Crystalline Silicon. Chem. Phys. Lett. 2001, 341, 185−194. (51) Murrell, J. N.; Laidler, K. J. Symmetries of Activated Complexes. Trans. Faraday Soc. 1968, 64, 371−377. (52) Nocedal, J. Updating Quasi-Newton Matrices with Limited Storage. Math. Comput. 1980, 35, 773−782. (53) Liu, D.; Nocedal, J. On the Limited Memory BFGS Method for Large Scale Optimization. Math. Prog. 1989, 45, 503−528. (54) Dijkstra, E. W. Numerische Math. 1959, 1, 269−271. (55) Prada-Gracia, D.; Gómez-Gardenes, J.; Echenique, P.; Fernando, F. Exploring the Free Energy Landscape: From Dynamics to Networks and Back. PLoS Comput. Biol. 2009, 5, 1−9. (56) Noé, F.; Krachtus, D.; Smith, J. C.; Fischer, S. Transition Networks for the Comprehensive Characterisation of Complex Conformational Change in Proteins. J. Chem. Theory Comput. 2006, 2, 840−857.

(57) Noé, F.; Fischer, S. Transition Networks for Modeling the Kinetics of Conformational Change in Macromolecules. Curr. Opin. Struct. Biol. 2008, 18, 154−162. (58) Becker, O. M.; Karplus, M. The Topology of Multidimensional Potential Energy Surfaces: Theory and Application to Peptide Structure and Kinetics. J. Chem. Phys. 1997, 106, 1495−1517. (59) Wales, D. J.; Miller, M. A.; Walsh, T. R. Archetypal Energy Landscapes. Nature 1998, 394, 758−760. (60) Evans, D. A.; Wales, D. J. Folding of the GB1 Hairpin Peptide from Discrete Path Sampling. J. Chem. Phys. 2004, 121, 1080−1090. (61) Wales, D. J. Calculating Rate Constants and Committor Probabilities for Transition Networks by Graph Transformation. J. Chem. Phys. 2009, 130, 204111-1−204111-7. (62) Wales, D. J. PATHSAMPLE: A Program for Generating Connected Stationary Point Databases and Extracting Global Kinetics. http://wwwwales.ch.cam.ac.uk/PATHSAMPLE/. (63) Marx, D.; Tuckerlan, M. E.; Hutter, J.; Parrinello, M. The Nature of the Hydrated Excess Proton in Water. Nature 1999, 397, 601−604. (64) Tuckerman, M.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab initio Molecular Dynamics Simulation of the Solvation and Transport of Hydronium and Hydroxyl Ions in Water. J. Chem. Phys. 1995, 103, 150−161. (65) CHEBI Database Entry:Creatine Zwitterion, http://www.ebi.ac. uk/chebi/searchId.do?chebiId=CHEBI:57947, 2012. (66) de Matos, P.; Alcántara, R.; Dekker, A.; Ennis, M.; Hastings, J.; Haug, K.; Spiteri, I.; Turner, S.; Steinbeck, C. Chemical Entities of Biological Interest: An Update. Nucleic Acids Res. 2010, 38, D249− D254. (67) Degtyarenko, K.; Hastings, J.; de Matos, P.; Ennis, M. ChEBI: An Open Bioinformatics and Cheminformatics Resource, Current Protocols in Bioinformatics; John Wiley and Sons, Inc.: New York, 2009, 26:14.9.1−14.9.20. (68) Degtyarenko, K.; de Matos, P.; Ennis, M.; Hastings, J.; Zbinden, M.; McNaught, A.; Alcántara, R.; Darsow, M.; Guedj, M.; Ashburner, M. ChEBI: a Database and Ontology for Chemical Entities of Biological Interest. Nucleic Acids Res. 2008, 36, D344−D350. (69) Wang, X.; Yin, Q. Determination of Isoelectric Point of Creatine by Ion Balance Method. J. Chem. Eng. Chin. Univ. 2003, 17, 569−574. (70) Lauzon, C. B.; van Zijl, P.; Stivers, J. T. Using the Water Signal To Detect Invisible Exchanging Protons in the Catalytic Triad of a Serine Protease. J. Biomol. NMR 2011, 50, 299−314. (71) Bai, Y.; Milne, J. S.; Mayne, L.; Englander, S. W. Primary Structure Effects on Peptide Group Hydrogen Exchange. Biochemistry 1972, 17, 75−86. (72) Zhou, J.; van Zijl, P. C. M. Chemical Exchange Saturation Transfer Imaging and Spectroscopy. Prog. Nucl. Magn. Reson. Spectrosc. 2006, 48, 109−136. (73) Sklenar, V. Gradient-Tailored Water Suppression for 1H-15N HSQC Experiments Optimized To Retain Full Sensitivity. J. Magn. Reson., Ser. A 1993, 102, 241−245. (74) Mori, S.; Abeygunawardana, C.; Berg, J. M.; van Zijl, P. C. M. NMR Study of Rapidly Exchanging Backbone Amide Protons in Staphylococcal Nuclease and the Correlation with Structural and Dynamic Properties. J. Am. Chem. Soc. 1997, 119, 6844−6852. (75) Barnett, C. B.; Wilkinson, K. A.; Naidoo, K. J. Molecular Details from Computational Reaction Dynamics for the Cellobiohydrolase I Glycosylation Reaction. J. Am. Chem. Soc. 2011, 133, 19474−19482. (76) Cui, Q.; Elstner, M.; Kaxiras, E.; Frauenheim, T.; Karplus, M. A QM/MM Implementation of the Self-Consistent Charge Density Functional Tight Binding (SCC-DFTB) Method. J. Phys. Chem. B 2001, 105, 569−585. (77) Meuwly, M.; Karplus, M. Simulation of Proton Transfer Along Ammonia Wires: an “ab initio” and Semiempirical Density Functional Comparison of Potentials and Classical Molecular Dynamics. J. Chem. Phys. 2002, 116, 2572−2585.

1975

dx.doi.org/10.1021/jp410172k | J. Phys. Chem. B 2014, 118, 1969−1975