Drastic Compensation of Electronic and Solvation Effects on ATP

Feb 21, 2017 - Hydrolysis of adenosine triphosphate (ATP) is the “energy source” for a variety of biochemical processes. In the present work, we a...
0 downloads 14 Views 2MB Size
Subscriber access provided by University of Newcastle, Australia

Article

Drastic Compensation of the Electronic and Solvation Effects on the ATP Hydrolysis Revealed Through a Large-Scale QM/MM Simulations Combined with a Theory of Solutions Hideaki Takahashi, Satoru Umino, Yuji Miki, Ryosuke Ishizuka, Shu Maeda, Akihiro Morita, Makoto Suzuki, and Nobuyuki Matubayasi J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b00637 • Publication Date (Web): 21 Feb 2017 Downloaded from http://pubs.acs.org on March 6, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Drastic Compensation of the Electronic and Solvation Effects on the ATP Hydrolysis Revealed through a Large-Scale QM/MM Simulations Combined with a Theory of Solutions

Hideaki Takahashi1*, Satoru Umino1, Yuji Miki1, Ryosuke Ishizuka2, Shu Maeda2, Akihiro Morita1, Makoto Suzuki3, and Nobuyuki Matubayasi2

1

Department of Chemistry, Graduate School of Science,

Tohoku University, Sendai Miyagi, 980-8578, Japan 2

Department of Chemical Engineering, Graduate School of

Engineering Science, Osaka University, Toyonaka Osaka, 560-8531, Japan 3

Department of Materials Processing, Graduate School of

Engineering, Tohoku University, Sendai Miyagi, 980-8579, Japan

* To whom correspondence should be addressed. electronic mail: [email protected], phone: +81-22-795-7722

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

Hydrolysis of adenosine triphosphate (ATP) is the “energy source” for a variety of biochemical processes. In the present work, we address key features of ATP hydrolysis: the relatively moderate value (about −10 kcal/mol) of the standard free energy ΔGhyd of reaction and the insensitivity of ΔGhyd to the number of excess electrons on ATP. We conducted quantum mechanical/molecular mechanical (QM/MM) simulation combined with the energy-representation theory of solutions to analyze the electronic-state and solvation contributions to ΔGhyd. It was revealed that the electronic-state contribution in ΔGhyd is largely negative (favorable) upon hydrolysis, due to the reduction of electrostatic repulsion accompanying the breakage of P-O bond. In contrast, the solvation effect was found to be strongly more favorable on the reactant side. Thus, we showed that a drastic compensation of the two opposite effects takes place, leading to the modest value of ΔGhyd at each number of excess electrons examined. The computational analyses were also conducted for pyrophosphate ions (PPi) and the parallelism between the ATP and PPi hydrolyses was confirmed. Classical molecular dynamics (MD) simulation was further carried out to discuss the effect of solvent environment; the insensitivity of ΔGhyd to the number of excess electrons was seen to hold in solvent water and ethanol.

2

ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. Introduction

Hydrolysis of adenosine triphosphate (ATP) is a vital chemical event to living systems 1-4

as a source of the free energy consumed by proteins in exhibiting their functions.

There

have been a lot of experimental and theoretical works to clarify the energetics of hydrolysis. kcal/mol

5-17

4,5

However, the underlying mechanism for the free energy release of about 10

associated with the P−O bond dissociation in ATP is not fully understood on the 1-3

molecular basis. In standard textbooks

it has been often speculated that the free energy

comes from the stabilization due to the fission of the highly repulsive phosphate ions and the electronic resonance effects in the resultant phosphate ions. However, the electronic and solvation free energies of the related ions, nor the possible competition between them have not been discussed at least on the quantitative basis. Thus, the energetics for such the fundamental reaction still remains unveiled. This is, of course, mainly due to the difficulty in the calculation of the free energy change associated with a chemical reaction of a large molecule in a condensed phase. Actually, it is only recently that the quantum mechanical 17

calculations were carried out for ATP

and its model systems

10,11,15,16

under the influence of

an explicit solvent represented with a classical force field. In this paper, we address the issue of the energetics of ATP hydrolysis by performing the large-scale quantum mechanical / molecular mechanical (QM/MM) simulations

18-23

where the whole ATP molecule is treated

as a quantum mechanical object. 5

Many years ago P. George and co-workers revealed experimentally that the free energy release is almost constant for the variation of the number Ne of the excess electrons on ATP. Explicitly, hydrolyses of ATP molecules with Ne = 3 and 4 yield free energies −11.4 and −10.8 kcal/mol, respectively. Such the constancy in the free energy release would be crucial for proteins to ensure their stable operations againt the change of pH in the environment. Our particular interest, therefore, will be placed on the constancy of the hydrolysis free energy insensitive to the charge state of ATP. To this end we analyze the free energy of hydrolyses of ATP and pyrophosphate ions often regarded as prototypes of ATP by decomposing the free energy into the contributions from the electronic-state and solvation effects. We also study the effect of the change in the solvent environment on the free energy of hydrolysis by performing a set of purely classical simulations.

24

2. Method For the evaluation of the free energy of a reaction in solution, a thermodynamic cycle which involves a reaction in the gas phase is often used. However, it is likely that the solutes with large Ne will not be stable in the absence of polar solvent. Thus, the calculated electronic energies of these solutes in the gas phase might be ambiguous. To circumvent this problem we directly evaluate the free energy ΔGhyd for the hydrolysis, e.g. ATP4− + H2O → ADP3− +

3

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

Pi− through the relation,

(

)

( )

(

)

ΔGhyd = G ADP 3− + G Pi− − G ATP 4− − G ( H 2O )

(1)

where the notation G(S) expresses the free energy of the solute S (= ATP4−, H2O, ADP3−, or Pi−) in water solution including the free energy of the electronic state. Then, each free energy G(S) can be decomposed into the contributions due to electronic and solvation free energies, thus, G ( S) = Gele ( S) + Gsol ( S)

(2)

where Gsol(S) denotes the solvation (hydration) free energy of the solute S with a fixed electron density

. The free energy due to the electronic energy including the electron

density fluctuation in solution constitutes the term Gele(S) in Eq. (2). Accordingly, ΔGhyd for the hydrolysis can be expressed in terms of the changes in these individual free energies, ΔGhyd = ΔGele + ΔGsol

.

(3)

Explicitly, ΔGsol for instance is expressed as

(

)

( )

(

)

ΔGsol = Gsol ADP 3− + Gsol Pi− − Gsol ATP 4− − Gsol ( H 2O )

.

(4)

The free energies Gsol(S) and Gele(S) are evaluated using a method, referred to as QM/MM-ER,

25,26

which

combines

the

energy-representation (ER) theory of solutions.

QM/MM

29-32

simulation

27,28

with

the

ER is an approximate scheme to compute

the solvation free energy using a functional constituting of distribution functions of solute-solvent interaction energy which is sampled during the course of simulation. Unlike standard free-energy methods such as free-energy perturbation, which introduce a number of intermediate states connecting the initial and final states of solute insertion, ER can provide the solvation free energy only with the simulations at the endpoints by formulating an approximate functional. The computational load is thus reduced, and the accuracy was assessed for the solvation free energies of a variety of solutes(15 amino-acid analogues), both hydrophobic and hydrophilic by making comparisons with the results given by numerically 33

exact approaches such as Bennet acceptance ratio (BAR) method.

It was revealed that the

difference between ER method and BAR employing the same force field is only ~1 kcal/mol for neutral species. It was, thus, suggested that the ER method is almost comparable in accuracy to other standard but numerically demanding approaches. Any form of solute-solvent potential can be adopted in this scheme, and in particular, the diffuse (cloud-like) nature of charge distributions in the QM region can be treated without any reduction, for example, to a set of point charges. In the framework of QM/MM-ER, the solvation free energy G(S) of a QM solute S is given by a sum of two contributions. The first

4

ACS Paragon Plus Environment

Page 5 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

is the free energy

due to two-body interaction between solute and solvent and the rest is

the many-body term δµ. In the present case these two contributions

and δµ are

equivalent to Gsol(S) and Gele(S) in Eq. (2), respectively. In the following we provide an outline for the QM/MM-ER method. We refer the readers to previous papers reviews

20-23

34-41

or

for more details.

The total energy Etot of a QM/MM system can be in general written as Etot = EQM [ X ] + EQM/MM [ X ] + EMM [ X ]

.

(5)

X in Eq. (5) stands for the set of the full coordinates for the solvent molecules, and EQM and EMM describe the energies of the QM and MM subsystems, respectively. EQM/MM is the interaction energy between QM and MM subsystems. In the QM/MM-ER approach, we configure the interaction EQM/MM ′ [ X ] as the sum of pairwise interactions, thus, (6)

where

represents the pair potential between the solute with a fixed density

and the ith solvent molecule with coordinate xi. Correspondingly, EQM ′ [ X ] is defined as the sum of EQM in Eq. (5) and the polarization energy due to the shift of the electron density from to n that depends on X, thus, .

In the present simulations the density

(7)

is set at the average electron density of each solute

in solution. The distributions ρ(ε) and ρ0(ε) of the pair interaction ε, which are constructed, respectively, in solution and pure solvent systems, serve as fundamental variables in evaluating the free energy Gsol(S) on the basis of the density functional theory of 29-32,42

solutions.

The energy η = EQM ′ [ X ] which involves the many-body interaction between

the solute and the solvent is responsible for the free energy Gele(S). The energy distribution functions P(η) and P0(η) are generated in solution and reference systems, respectively. Then, the free energy functional in Ref. 26 is applied to these distributions to yield Gele(S). The distribution functions for the free energies Gele(S) and Gsol(S) are also constructed through QM/MM simulations. 3. Computational Details We perform a set of QM/MM (quantum mechanical / molecular mechanical) simulations which employ the real-space grids (RSG)

43-47

to achieve a high parallel efficiency in the

Kohn-Sham density functional theory (KS-DFT)

48

calculations for the QM region. The 20,23

methodological details for KS-DFT using RSG are presented in our reviews.

5

ACS Paragon Plus Environment

An ATP

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

molecule, regarded as a QM solute, is embedded in a spherical water droplet with radius of 25 Å, where the mass center of the ATP molecule is placed at the center of the sphere. The water droplet contains 2160 molecules represented with the SPC/E model.

49

The water molecules

near the surface of the sphere are repelled by a LJ potential. The ATP molecules with Ne = 0, 3 and 4 are treated in the present simulations. The molecular geometries of these solutes are optimized by conducting the Gaussian 09 (G09) package B3LYP

51-53

/aug-cc-pVDZ

54

50

with the theoretical level of

We employ the polarizable continuum model (PCM)

55

to realize

the influence of the water solvent. The optimized structures for the related species are illustrated in Table S1 in Supporting Information. The specific hydrogen bond interactions between solute and solvent cannot be realized in the PCM approach, the stable structures in the PCM might be different from those in the actual bulk water. It is naturally expected, however, that destabilization in the solute-solvent interaction can be compensated with the stabilization in the solvent-solvent interaction in many cases. Thus, the errors in the free energies which arises from the PCM geometries will not be so serious. A QM/MM snapshot is presented in Fig. 1 to show graphically the relative size of the droplet to ATP. ----------------------------------------Figure 1 ----------------------------------------During the QM/MM simulations the geometry of the solute is kept fixed and its electronic state is described with KS-DFT. We also take into consideration the free energies due to the rotational and vibrational motion of each solute. We evaluate these thermal contributions as well as the zero-point energies by invoking Freq option in G09 package. The solute wave functions are enclosed within a rectangular QM cell of which x, y, and z dimensions are 26.7 Å, 18.2 Å, and 12.1 Å, respectively. The grid points are uniformly placed along each axis 51,53

with the width of 0.152 Å. The BLYP functional

is used to calculate the exchange

correlation energy of the solute. All the QM/MM simulations for ATP are performed on the massively parallel computers(Cray XC30, SX-ACE). For the parallel KS-DFT calculations the QM cell is typically divided into 128 rectangular sub-cells by slicing the cell into 8, 4 and 4 layers along the x, y, and z directions, respectively. Then, the computational task associated with each sub-cell is assigned to a node in a parallel computer. The data communication among these sub-cells are commanded using the MPI libraries. The specific ‘do-loops’ which consume a lot of CPU time in the code are also parallelized with 4 cores on each node by conducting the OpenMP directives. Consequently, 1 SCF (self-consistent field) step for KS-DFT can be processed within 0.6 sec. The average electron distribution n! for each solute is constructed through a 50-ps QM/MM simulation. For the constructions of the distribution functions ρ(ε) and ρ0(ε) for the free energy Gsol(S) we devote typically 200-ps and 400-ps QM/MM simulations, respectively, while 100-ps simulations are used to yield the distribution functions for the electronic free energy

6

ACS Paragon Plus Environment

Page 7 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Gele(S). To ensure the sufficient overlap between the distributions ρ(ε) and ρ0(ε) we perform 1.7, 0.8, and 4.8 ns simulations to build the distribution functions ρ0(ε) for the solutes ATP4−, ATP3−, and ADP3−, respectively. We construct the distribution functions for the solvent molecules within the sphere Ω of radius a0 = 23 Å of which center is set at the mass center of 56

the solute. For ionic solutes we incorporate the Born correction

to Gsol(S) to evaluate the

electrostatic contribution from the water molecules outside the sphere Ω. Explicitly, the Born equation

56

is given by −165.9×Ne2(1-1/ε)/a0 in the unit of kcal/mol where a0 is the radius of

the ion cavity and ε is the dielectric constant. The adequacy of the Born estimation is fully discussed in Supporting Information. Since the pyrophosphate ion (PPi) has often been considered as a model system of ATP, another part of our calculation is devoted to the QM/MM simulations for the hydrolyses of the PPi (Ne = 0, 1, 2, 3, and 4) in water. The geometries of these solutes are optimized in the same way as those of APT molecules. For these relatively small solutes we employed a cubic MM simulation cells with periodic boundary conditions. For PPi with charges 0 and −1 we employ MM cell of the size L = 24.6 Å which contains 494 SPC/E water molecules, while we employ larger cell (L = 49.3 Å) containing 3994 solvent molecules for the systems of Ne = 2, 3, and 4. The size of the cubic QM cell for pyrophosphate is set at 15 Å. The grid width for real-space cell is the same as that used for the ATP simulations. The distribution functions are constructed for the solvent molecules within the sphere Ω' of radius a0' = 22 Å of which center is set at the mass center of the solute. For the ionic solutes in the PPi reactions, the free energy contribution from the water molecules outside the sphere Ω' is also evaluated by Born equation. The LJ parameters in the OPLS force field

57

are assigned to the interaction sites in ATP

and PPi. Then, some of the size parameters on oxygen atom are modified so that the QM/MM pair potential between solute and solvent reproduces the corresponding QM potential. The parameters employed in the present simulations and the details of the modification are provided in Table S1. The Newtonian equations of motion for the MD simulation are numerically solved using the velocity Verlet algorithm

24

where the time step is set at 1 fs and

the internal coordinates of the solvent molecules are fixed. 1-fs time step would be adequate to describe translational and rotational motions of the solvent molecules. Throughout the present simulations the thermodynamic condition is set at ρ = 1.0 g/cm3 and T = 300 K. The reactions of PPi as well as ATP treated in this work are summarized in Table 1 for later references. The other miscellaneous details for the QM/MM-ER simulations are provided also in Supporting Information. ----------------------------------------Table 1 -----------------------------------------

7

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

4. Results and Discussion The distribution ρ(ε) of the pair interaction ε between ATP4− and solvent molecules in the solution system are presented in Fig. 2(a). Also shown is the energy distribution ρ0(ε) in the pure solvent system where the solvent configurations are generated under the condition that solute-solvent interaction is absent. ----------------------------------------Figures 2(a) and 2(b) ----------------------------------------We observe in Fig. 2(a) a distinct population of the solute-solvent pair interaction over the energy coordinate of −35 ~ −15 kcal/mol, though not conspicuous. The average amount of ESP charges on the tri-phosphate moiety ((PO3)3-O-) of ATP4− in aqueous solution is evaluated as −4.23 e, which indicates that the excess electrons are confined within the relatively small atomic group irrespective of the large repulsive forces among the electrons. The destabilization associated with the electron repulsion may be compensated by hydration effect. The inset in Fig. 2(a) shows the superimposed snapshots of the water molecules of which solute-solvent interaction ε are within the range of −35 < ε < −30 kcal/mol and within −30 < ε < −25 kcal/mol. As expected, the water molecules with strong binding energies for the solute are localized around the tri-phosphate moiety. Here, we note that the water molecules with the interaction energy in the range of −35 < ε < −30 kcal/mol form typically two hydrogen bondings (HB) with the oxygen atoms of the Pi groups, while those in the range of −30 < ε < −25 kcal/mol are mostly connected with single HBs with the oxygen atoms in the phosphate groups. We apply an approximate free energy functional developed in Ref. 25, 30-32 to these energy distributions and the correlation matrix in the pure solvent system, which yields the solvation free energy Gsol(ATP4−) = −658.2 kcal/mol including the Born correction (−113.9 kcal/mol). Figure 2(b) shows the distribution functions for Gsol(ADP3−) in the solution and the pure solvent systems. It is recognized that the energy distributions for ADP3− span over narrower ranges of ε as compared with those for ATP4− as a result that the solute-solvent interaction is decreased. Specifically, the population around −35 kcal/mol is decreased in the solution system of ADP3− as an implication that the solute-solvent interaction is weakened. With these distributions we obtain the solvation free energy Gsol(ADP3−) = −430.0 kcal/mol that includes the Born correction (−64.1 kcal/mol). To calculate the electronic free energy Gele(ATP4−) the distributions of the energy η = EQM ′ [ X ] defined in Eq. (7) are generated in the solution and in the reference systems. We note that η includes the distortion energy of the electronic energy of the QM solute in solution. In the solution system the solute fully couples with solvent molecules, while in the reference system the solute electron density is being fixed at

. The distributions P(η) in the solution

and P0(η) in the reference system for ATP4− are shown in Fig. 3(a), where the average = −339.94113 a.u. in the solution is taken as the standard for the η value.

8

ACS Paragon Plus Environment

Page 8 of 29

Page 9 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

----------------------------------------Figures 3(a) and 3(b) ----------------------------------------It is observed in the figure that the peak of P(η) in the solution system shifted toward the lower energy region as compared to P0(η) implying the stabilization of the system due to the electron density fluctuation of the solute around

. By substituting these functions P(η) and

P0(η) to the free energy functional defined in Eq. in Ref. 26 we evaluate the electronic free energy as Gele(ATP4−) = −15.8 kcal/mol when expressed relative to the average = −339.94113 a.u. for ATP4− in solution. We note that includes the thermal free energy and the zero-point vibrational correction. Thus, the total free energy G(ATP4−) is given as −674.0 kcal/mol by the sum of the free energies Gsol(ATP4−) and Gele(ATP4−), where the energy is similarly taken as the standard. The distributions P(η) and P0(η) of the energy η = EQM ′ [ X ] for the solute ADP3− are presented in Fig. 3(b). The overall behaviors of these functions are common to those of the solute ATP4−. However, the peak positions of P(η) and P0(η) of ADP3− are shifted toward larger energy coordinate by 2~3 kcal/mol. This can be attributed to the fact that ADP3− is less polarizable as compared to ATP4−. Thus, the stabilization due to the density polarization of ADP3− in solution is relatively small as compared to ATP4−. As a result the electronic and total free energies Gele(ADP3−) and G(ADP3−) are evaluated as −13.5 kcal/mol and −443.5 kcal/mol, respectively, relative to = −285.69306 a.u. With the similar manner we computed the free energies of the other species (Pi1− and H2O) relevant to the hydrolysis of ATP4−. The results are summarized in Table 2. ----------------------------------------Table 2 ----------------------------------------We note that the value in the parenthesis in Gele shows the zero-point corrected for each solute and is used as the standard of the free energy Gele. Substituting these values to Eq. (1) we obtain the free energy ΔGhyd of hydrolysis of ATP4− as −16.3 kcal/mol showing a rather good agreement with the experimental value of −10.7 kcal/mol. Decomposition of

ΔGhyd into the electronic and solvation free energies, ΔGele and ΔGsol, as expressed in Eq. (3) provides an intriguing mechanism for the ATP hydrolysis. ΔGele is evaluated as −172.5 kcal/mol and, thus, the electronic contribution is found to be highly exothermic, while the change in the solvation free energy is endothermic (ΔGsol = 156.1 kcal/mol) on the contrary. It is, thus, revealed that a drastic cancellation of these opposite effects results in the modest free energy release ΔGhyd of ATP in aqueous solution. The exothermicity in ΔGele can be merely attributed to the decrease in the electronic repulsion as a result of the fragmentation of the highly negative solute. The endothermicity of ΔGsol can also be understood on the qualitative basis by consulting the Born equation which describes the electrostatic solvation free energy

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of an ion. As shown in the Born equation the electrostatic free energy is proportional to the square of the charge, and the fragmentation of an ion into pieces with smaller absolute charges will cause substantial destabilization in the solvation free energy. In addition, the total volume of the ionic fragments will be larger than that of the original solute, which also causes the destabilization in ΔGsol. ----------------------------------------Figure 4 ----------------------------------------To examine the cancellation effect in other systems with different excess electrons we plot in Fig. 4 the free energy changes ΔGhyd, ΔGele and ΔGsol with respect to the variation of Ne for pyrophosphate ions (PPi) as well as ATPs. It is quite striking in the curves for PPi that the free energies ΔGele largely compensate for the destabilized ΔGsol resulting in modest ΔGhyd irrespective of the value of Ne. Interestingly, at Ne = 1 the effects of ΔGele and ΔGsol are inverted. As a specific feature in the case of Ne = 1 the cavity volume of the ion is decreased after fragmentation since one of the resulting fragments is neutral, which causes the stabilization in the free energy change ΔGsol and simultaneously the destabilization in ΔGele. Anyway, our computational results for PPi hydrolyses reasonably realize the constancy of the free energy release ΔGhyd against the value of Ne. Explicitly, ΔGhyd for Ne = 0 ~ 4 lie within the range from −7.3 to −17.8 kcal/mol. In Table 3 the values of ΔGhyd given by the QM/MM-ER method are compared with the experimental data ranging from −7.1 to −10.4 kcal/mol. We also examined ΔGsol due to solvation effect by using the Born equation for which radius a0 of each ion in PPi reactions is determined by the ‘volume’ option in the Gaussian 09 package under the PCM environment. The results from the Born estimate are plotted in Fig. 4. Regardless of its rude approximation, it is found that the Born equation reproduces the proper behavior of ΔGsol at least qualitatively, which substantiates our interpretation given above. Specifically, at Ne = 3 and 4 the Born approximation favorably reproduces the results of ΔGsol given by QM/MM-ER suggesting that the long range Coulomb interactions for the solutes with larger Ne dominate the solvation free energies. We also make a remark on the overestimations of ΔGhyd observed as a common tendency in the hydrolyses of ATP and PPi. In this work we employed the GGA functional (BLYP) in the KS-DFT calculation for the QM solutes. It is well known that a conventional GGA functional suffers 58

from the delocalization error

for excess charges on anions in particular. The erroneously

delocalized charges on the solutes will lead to the destabilization in the solute-solvent interactions. This effect is larger in the reactant states with larger excess charges and hence stabilizes relatively the product states. Thus, we consider the major source of errors in ΔGhyd comes from delocalization error in the GGA functional in addition to the inadequacy of the force field.

10

ACS Paragon Plus Environment

Page 10 of 29

Page 11 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

----------------------------------------Table 3 ----------------------------------------The free energy changes ΔGele and ΔGsol for the reactions of ATP molecules with Ne = 0, 3, and 4 show the similar trend with those for PPi. However, absolutes of them are smaller than those of the corresponding PPi reactions. In the ATP reactions the excess electrons (Ne = 3 and 4) are on three phosphate groups, while they are confined within two groups in the PPi reactions. Consequently, the solvation effects as well as the electron repulsions are decreased, which may also lead to the reduction of the absolutes of the free energy changes ΔGele and

ΔGsol associated with the hydrolyses. We, thus, find that the constancies in ΔGhyd for PPi and ATP hydrolyses are the consequences of the drastic cancellations between free energy changes ΔGele and ΔGsol associated with chemical bond rearrangements. Here, we remind that it is often mentioned in the textbooks that the major sources of the free energy associated with the P−O bond dissociation are (i) the relaxation of the electronic repulsion and (ii) the stabilization of the solvation free energies of the relevant species after the hydrolyses. As a major conclusion of this paper we stress that the solvation effect actually destabilizes the product state, as completely opposed to the notion (ii). In fact the solvation effect almost cancels the large stabilization in the electronic contribution resulting in the modest free energy release. Our final issue is to discuss the dependence of ΔGsol on the solvent environment. To this end, we performed additional sets of purely classical MD simulations and evaluated the solvation free energies Gsol of the reactive species in solvent water and in solvent ethanol. In these calculations, the solute structures were frozen and were set to be identical to those used in the QM/MM simulations. The Lennard-Jones parameters were also the same, and the partial charges on the atomic sites were determined with the ESP procedure for which details are described in Supporting Information. The solvation free energies were then computed with the energy-representation method,29-32 and all the MDs were conducted with GROMACS 5.0.5 in the NVT ensemble at 300 K with periodic boundary condition and minimum image convention.59,60 The number of solvent molecules in the MD unit cell was 5000 when the solvent is water, and was 1500 when the solvent is ethanol. The water and ethanol molecules were modeled with SPC/E49 and OPLS-AA57, respectively, and the cell size was set so that the solvent density is equal to 0.997 and 0.785 g/cm3 for water and ethanol, respectively. All the MD simulations were carried out for 10 ns with a sampling interval of 0.1 ps. The leapfrog stochastic dynamics integrator was employed at a time step of 2 fs and an inverse friction constant of 0.2 ps.61 The electrostatic interaction was treated by the smooth particle-mesh Ewald (PME) method.62 The PME parameters were 12 Å for the real-space cutoff, 4 for the spline order, 10-5 for the relative tolerance, and 53 for the reciprocal-space mesh size along each of the x, y, and z directions. The Lennard-Jones (LJ) interactions were

11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

spherically truncated at 12 Å, with switching range of 10−12 Å.63 Both of the real-space part of PME interaction and the LJ interaction were truncated on the atom-atom basis, and the LJ interaction was not subject to the long-range correction. The self energy of PME interaction, which refers to the interaction of the solute with its own images and neutralizing background, was added to the solvation free energy to correct the finite-size effect.64 The lengths of all the bonds with hydrogen atoms in ethanol were fixed with LINCS,65 and the SETTLE algorithm was used to keep the water molecules rigid.66 ----------------------------------------Figure 5 ----------------------------------------The solvation free energies thus computed with classical MD are listed in Table 2 and Table S2 in the Supporting Information, and in Fig. 5 we show the changes in the solvation free energies ΔGsol in the three reactions of R1, R2, and R3 in Table 1; their numbers of excess electrons Ne are 0, 3, and 4, respectively. The electronic state of the solute does not fluctuate in the classical MD, and the computed results for solvation in Table 2 and in Fig. 5 are to be compared to the corresponding, two-body terms Gsol in the QM/MM calculations. It is seen in Table 2 and Table S2 in the Supporting Information that when the solvent is fixed to water, the solvation free energy obtained from classical MD is larger in magnitude than that 28,35

from QM/MM. This is in line with previous observations for anions,

and is due to the

reduction of diffuse electronic distribution into a set of point charges. The discrepancy between the classical and QM/MM results are much smaller for ΔGsol in Fig. 5, though, and when the electronic part is taken from QM/MM and is combined without modification to the classical results, the free energy change ΔGhyd of hydrolysis is found to be −19.8, −15.6, and 2.9 kcal/mol at Ne = 0, 3, and 4, respectively. The extent of constancy of the hydrolysis free energy ΔGhyd is thus comparable between QM/MM and classical MD in solvent water. The solvation free energies in solvent ethanol are also listed in Table 2. They are smaller in magnitude than in solvent water, showing the sensitivity of the extents of stabilization especially for highly charged species. Still, Fig. 5 shows that the difference between solvent water and ethanol is much reduced for the solvent effect in the hydrolysis free energy ΔGsol. When the electronic part computed in water is taken from QM/MM without modification, as done in the cases of classical MD in solvent water, ΔGhyd in solvent ethanol remains at −15.8, −28.9, and −40.5 kcal/mol with Ne = 0, 3, and 4, respectively. This extent of constancy is comparable to the water results described above. The large values of the solvation free energies in Table 2 are due to electrostatic interaction between the solutes and the polar solvents, and those large effects are well cancelled for the excess-electron dependence of the ATP hydrolysis. The remarkable insensitivity of ΔGhyd to the excess electron is thus shown in a solvent environment with reduced dielectric constant. It should be noted, however, that

ΔGhyd would be highly exothermic in a condition where dielectric constant of the environment

12

ACS Paragon Plus Environment

Page 13 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

is much smaller than that of the solvent water. 5. Conclusion We performed large scale QM/MM-ER simulations for the hydrolyses of ATP as well as pyrophosphate ions in aqueous solutions to clarify the mechanism for the constancy of the free energy release ΔGhyd associated with the reactions irrespective of the number Ne of the excess electrons on the solutes. For this purpose we decomposed the free energy change

ΔGhyd into the contributions ΔGele and ΔGsol due, respectively, to the electronic energy and the solvation effect by means of the theory of solutions. We found that the solvation part ΔGsol substantially destabilizes as Ne increases, while ΔGele stabilizes in contrast. The constancy of the free energy ΔGhyd for the variation of Ne is realized by the fact that ΔGele almost compensates for the destabilized ΔGsol resulting in modest ΔGhyd. We also found the same trend in the energetics of the hydrolyses of pyrophosphoric acids often regarded as the prototype of ATP. The present results show clear contrast to the description in the standard text books which speculate that the free energy release of about 10 kcal/mol associated with the P-O bond dissociation in ATP is partly due to the stabilization in ΔGsol. We also examined the effect of the solvent environment by performing additionally the classical MD simulations for solvent ethanol. It was revealed that the free energy gains for the ATP hydrolyses are rather insensitive to the change in the dielectric constant of the solvent, showing the stability of the free energy release against the variation of the dielectric constant of the solvent. Supporting Information: ・Table S1: LJ parameters used in the QM/MM simulations and ESP charges employed in the full classical simulations. ・Table S2: Free energy components for the ATP hydrolyses. ・Table S3: Free energy components for the PPi hydrolyses. ・Miscellaneous Computational Details ・Effect of Born Correction ・Statistical Errors of Free-Energy Calculations

ACKNOWLEDGEMENTS This work is supported by the Grants-in-Aid for Scientific Research on Innovative Areas (Nos. 20118002, 20118008, and 23118701) and the Elements Strategy Initiative for Catalysts and Batteries from the Ministry of Education, Culture, Sports, Science, and Technology, by the Grants-in-Aid for Scientific Research (No. 15K13550, 25620004, and 26240045) from the Japan Society for the Promotion of Science, and by Computational Materials Science Initiative, Theoretical and Computational Chemistry Initiative, and Strategic Programs for 13

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Innovative Research of the Next-Generation Supercomputing Project. The simulations were conducted partly using COMA at University of Tsukuba and Cray XC30 at Kyoto University through the HPCI System Research Project (Project IDs: hp160007, hp160013, hp160019, and hp160214).

References 1

Voet, D.; Voet, J. G.; Pratt, C. W. Fundamentals of Biochemistry 4th ed; Wiley: Hoboken, 2013.

2

Nelson, D. L.; Cox, M. M. Lehninger Principles of Biochemistry 6th ed; W. H. Freeman and Company: New York, 2013.

3

Berg, J. M.; Tymoczko, J. L.; Gatto, Jr. G. L.; Stryer, L. Biochemistry 8th ed; W. H. Freeman and Company: New York, 2015.

4

Meyerhof, O.; Lohmann, K. Über Energetische Wechselbeziehungen Zwischen dem Umsatz der Phosphorsäureester im Muskelextrakt. Biochem. Z. 1932, 253, 431-461.

5

George, P.; Witonsky, R. J.; Trachtman, M.; Wu, C.; Dorwart, W.; Richman, L.; Richman, W.; Shurayh, F.; Lentz, B. “Squiggle-H2O”. An Enquiry into the Importance of Solvation Effects in Phosphate Ester and Anhydride Reactions. Biochim. Biophys. Acta 1970, 223, 1−15.

6

Kodama, T. Thermodynamic Analysis of Muscle ATPase Mechanisms. Physiol. Rev. 1985, 65, 467-551.

7

Alberty, R. A.; Goldberg, R. N. Standard Thermodynamic Formation Properties for the Adenosine 5’-Triphosphate Series. Biochemistry 1992, 31, 10610−10615.

8

Pepi, F.; Ricci, A.; Rosi, M.; Di Stefano, M. Gaseous H5P2O8- Ions: A Theoretical and Experimental Study on the Hydrolysis and Synthesis of Diphosphate Ion. Chem. Eur. J. 2004, 10, 5706-5716.

9

Colvin, M. E.; Evleth, E.; Akacem, Y. Quantum Chemical Studies of Pyrophosphate Hydrolysis. J. Am. Chem. Soc. 1995, 117, 4357−4362.

10

Klähn, M.; Rosta, E.; Warshel, A. On the Mechanism of Hydrolysis of Phosphate Monoesters Dianions in Solutions and Proteins. J. Am. Chem. Soc. 2006, 128, 15310−15323.

11

Grigorenko, B. L.; Rogov, A. V.; Nemukhin, A. V. Mechanism of Triphosphate Hydrolysis in Aqueous Solution: QM/MM Simulations in Water Clusters. J. Phys. Chem. B 2006, 110, 4407−4412.

12

Ross, J. Energy Transfer from Adenosine Triphosphate. J. Phys. Chem. B 2006, 110, 6987− 6990.

13

Ruben, E. A.; Plumley, J. A.; Chapman, M. S.; Evanseck, J. D. Anomeric Effect in

14

ACS Paragon Plus Environment

Page 14 of 29

Page 15 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

“High Energy” Phosphate Bonds. Selective Destabilization of the Scissile Bond and Modulation of the Exothermicity of Hydrolysis. J. Am. Chem. Soc. 2008, 130, 3349−3358. 14

Arabi, A. A.; Matta, C. F. Where is Electronic Energy Stored in Adenosine Triphosphate? J. Phys. Chem. A 2009, 113, 3360−3368.

15

Kamerlin, S. C. L.; Warshel, A. On the Energetics of ATP Hydrolysis in Solution. J. Phys. Chem. B 2009, 113, 15692−15698.

16

Hong, J.; Yoshida N.; Chong, S.-H. Lee, C.; Ham, S.; Hirata, F. Elucidating the Molecular Origin of Hydrolysis Energy of Pyrophosphate in Water. J. Chem. Theory Comput. 2012, 8, 2239-2246.

17

Wang, C.; Huang, W.; Liao, J.-L. QM/MM Investigation of ATP Hydrolysis in Aqueous Solution. J. Phys. Chem. B 2015, 119, 3720− 3726.

18

Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; Wiley: New York, 1991.

19

Gao, J.; Xia, X. A Priori Evaluation of Aqueous Polarization Effects Through Monte Carlo QM-MM Simulations. Science 1992, 258, 631-635.

20

Combined QM/MM calculations in chemistry and biochemistry in J. Mol. Struct.: THEOCHEM Vol. 632 pp 1-336 (Special Issue); Ruiz-López M. F. Eds.; Elsevier: Amsterdam, 2003.

21

Solvation effects on Molecules and Biomolecules: Challenges and Advances in Computational Chemistry and Physics, Vol. 6; Canuto, S. Eds.; Springer: Heidelberg, 2008.

22

Combining Quantum Mechanics and Molecular Mechanics: Some Recent Progresses in QM/MM Methods, Advances in Quantum Chemistry Vol. 59; Canuto, S. Eds.; Elsevier: Amsterdam, 2010.

23

Quantum Modeling of Complex Molecular Systems, Challenges and Advances in Computational Chemistry and Physics, Vol. 21; Rivail, J.-L.; Ruiz-Lopez M.; Assfeld X. Eds.; Springer: Heidelberg, 2015.

24

Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987.

25

Takahashi, H.; Matubayasi, N.; Nakahara, M.; Nitta, T. A Quantum Chemical Approach to the Free Energy Calculations in Condensed Systems: The QM/MM Method Combined with the Theory of Energy Representation. J. Chem. Phys. 2004, 121, 3989-3999.

26

Takahashi, H.; Omi, A.; Morita, A.; Matubayasi, N. Simple and Exact Approach to the Electronic Polarization Effect on the Solvation Free Energy: Formulation for Quantum-Mechanical/ Molecular-Mechanical System and its Applications to Aqueous Solutions. J. Chem. Phys. 2012, 136, 214503.

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

27

Takahashi, H.; Hori, T.; Hashimoto, H.; Nitta, T. A Hybrid QM/MM Method Employing Real Space Grids for QM Water in the TIP4P Water Solvents. J. Comput. Chem. 2001, 22, 1252-1261.

28

Takahashi, H.; Takei, S.; Hori, T.; Nitta, T. A Real-Space-Grid QM/MM Study on the Ionic/Radical Association Reaction in Aqueous Phase: HCHO+OH→HCHO−OH J. Mol. Struct.: THEOCHEM, 2003, 632, 185-195.

29

Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energetic Representation. I. Formulation. J. Chem. Phys. 2000, 113, 6070-6081.

30

Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energy Representation. II. Functional for the Chemical Potential. J. Chem. Phys. 2002, 117, 3605-3616; Erratum. J. Chem. Phys. 2003, 118, 2446.

31

Matubayasi, N.; Nakahara, M. Theory of Solutions in the Energy Representation. III. Treatment of the Molecular Flexibility. J. Chem. Phys. 2003, 119, 9686-9702.

32

Sakuraba, S; Matubayasi, N. ERmod: Fast and Versatile Computation Software for Solvation Free Energy with Approximate Theory of Solutions. J. Comput. Chem. 2014, 35, 1592-1608.

33

Karino, Y; Fedorov, M. V; Matubayasi, N. End-point Calculation of Solvation Free Energy of Amino-acid Analogs by Molecular Theories of Solution. Chem. Phys. Lett. 2010, 496, 351-355.

34

Takahashi, H.; Satou, W.; Hori, T.; Nitta, T. An Application of the Novel Quantum Mechanical/Molecular Mechanical Method Combined with the Theory of Energy Representation: An Ionic Dissociation of a Water Molecule in the Supercritical Water. J. Chem. Phys. 2005, 122, 044504.

35

Takahashi, H.; Kawashima, Y.; Nitta, T.; Matubayasi, N. A Novel Quantum Mechanical / Molecular Mechanical Approach to the Free Energy Calculation for Isomerization of Glycine in Aqueous Solution. J. Chem. Phys. 2005, 123, 124504.

36

Hori, T.; Takahashi, H.; Nakano, M.; Nitta, T.; Yang, W. A QM/MM Study Combined with the Theory of Energy Representation: Solvation Free Energies for Anti / Syn Acetic Acids in Aqueous Solution. Chem. Phys. Lett. 2005, 419, 240-244.

37

Hori, T.; Takahashi, H.; Furukawa, S.; Nakano, M.; Yang, W. Computational Study on the Relative Acidity of Acetic Acid by the QM/MM Method Combined with the Theory of Energy Representation. J. Phys. Chem. B 2007, 111, 581-588.

38

Takahashi, H.; Ohno, H.; Kishi, R.; Nakano, M.; Matubayasi, N. Computation of the Reduction Free Energy of Coenzyme in Aqueous Solution by QM/MM-ER Method. Chem. Phys. Lett. 2008, 456, 176-180.

39

Takahashi, H.; Miki, F.; Ohno, H.; Kishi, R.; Ohta, S.; Furukawa, S.; Nakano, M. J. Math. Chem. Hydration Effects on the Reaction with an Open-Shell Transition State: QM/MM-ER Study for the Dehydration Reaction of Alcohol in Hot Water. 2009, 46,

16

ACS Paragon Plus Environment

Page 16 of 29

Page 17 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

781-794. 40

Takahashi, H.; Maruyama, K.; Karino, Y.; Morita, A.; Nakano, M.; Jungwirth, P.; Matubayasi, N. Energetic origin of Proton Affinity to the air/water interface. J. Chem. Phys. B 2011, 115, 4745-4751.

41

Takahashi, H.; Iwata, Y.; Kishi, R.; Nakano, M. The QM/MM-ER studies for the origin of the antioxidative properties of MCI-186 in aqueous solutions. Int. J. Quantum Chem. 2011, 111, 1748-1762.

42

Hansen, P.; McDonald, I. R. Theory of Simple Liquids, 3rd ed.; Academic: London, 2006.

43

Chelikowsky, J. R.; Troullier, N.; Saad, Y. Finite-Difference-Pseudopotential Method: Electronic Structure Calculations without a Basis. Phys. Rev. Lett. 1994, 72, 1240-1243.

44

Chelikowsky, J. R.; Troullier, N.; Wu, K.; Saad, Y. Higher-Order Finite-Difference Pseudopotential Method: An Application to Diatomic Molecules. Phys. Rev. B 1994, 50, 11355-11364.

45

Takahashi, H.; Hori, T.; Wakabayashi, T.; Nitta, T. Chem. Lett. A Density Functional Study for Hydrogen Bond Energy by Employing Real Space Grids. 2000, 29, 222-223.

46

Takahashi, H.; Hori, T.; Wakabayashi, T.; Nitta, T. Real Space Ab Initio Molecular Dynamics Simulations for the Reactions of OH Radical / OH Anion with Formaldehyde. J. Phys. Chem. A 2001, 105, 4351-4358.

47

Takahashi, H.; Hori, T.; Hashimoto, H.; Nitta, T. A Hybrid QM/MM Method Employing Real Space Grids for QM Water in the TIP4P Water Solvents. J. Comp. Chem. 2001, 22, 1252-1261.

48

Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation effects. Phys. Rev. 1965, 140, A1133-A1138.

49

Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269-6271.

50

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al, Gaussian 09, Revision C.01; Gaussian, Inc.: Wallingford, CT, 2010.

51

Becke, A. D. Density-functional Exchange-energy Approximation in Calculations with Correct Asymptotic-behavior. Phys. Rev. A 1988, 38, 3098-3100.

52

Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648-5652.

53

Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-energy Formula into a Functional of The Electron Density. Phys. Rev. B 1988, 37, 785-789.

54

Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007-1023.

55

Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation

17

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Models. Chem. Rev. 2005, 105, 2999-3094. 56

Born, M.; Volumen und hydratationswarme der ionen. Z. Phys. 1920, 1, 45-48.

57

Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of The OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236.

58

Cohen, A. J; Mori-Sanchez, P; Yang W. Insights into Current Limitations of Density Functional Theory. Science, 2008, 321, 792-794.

59

Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435-447.

60

Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High Performance Molecular Simulations through Multi-level Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1-2, 19-25.

61

van Gunsteren, W. F.; Berendsen, H. J. C. A Leap-frog Algorithm for Stochastic Dynamics.

62

Mol. Simul. 1988, 1, 173–185.

Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593.

63

van der Spoel, D.; van Maaren, P. J. The Origin of Layer Structure Artifacts in Simulation of Liquid water. J. Chem. Theory Comput. 2006, 2, 1-11.

64

Hummer, G.; Pratt, L. R.; García, A. E. Free Energy of Ionic Hydration. J. Phys. Chem. 1996, 100, 1206-1215.

65

Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. 3 LINCS: a Linear constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463−1472.

66

Miyamoto, S.; Kollman, P. A.: SETTLE: an Analytical Version of the SHAKE and RATTLE Algorithms for Rigid Water Models. J. Comput. Chem. 1992, 13, 952-962.

18

ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1.

Hydrolyses for ATP and PPi treated in the present QM/MM-ER

simulations. Hydrolyses of adenosine triphosphate (ATP): R1: R2: R3 Hydrolyses of pyrophosphate (PPi): R4 R5 R6 R7 R8

Table 2. The free energy components in the unit of kcal/mol for the species related to the hydrolysis reaction R3 presented in Table 1. The value in the parenthesis is the average of the energy EQM in Eq. (7) used as the standard of the free energy Gele for each solute. EQM is corrected with the zero-point energy and includes the thermal (rotational and vibrational) free energies. EQM are given in the atomic unit. The notation, e.g. Classical MD(ethanol) indicates that a purely classical molecular dynamics simulation is performed for the ethanol solvent.

Ne = −4

Gele kcal/mol (a.u.) QM/MM (water)

Gsol kcal/mol QM/MM (water)

Gsol kcal/mol Classical MD (water)

Gsol kcal/mol Classical MD (ethanol)

ATP4−

−15.8 (−339.94113)

−658.2

−838.4

−642.3

H 2O

−1.5 (−17.08730)

−8.6

−7.9

−6.2

ADP3−

−13.5 (−285.69306)

−430.0

−558.8

−431.6

Pi−

−4.5 (−71.60914)

−80.6

−112.0

−85.1

19

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 3.

Page 20 of 29

Comparison of the free energies ΔGhyd of the hydrolysis reactions

R1~R8 listed in Table1. The free energy values given by the present work (PW) can be constructed from the free energy components in the Tables S2 and S3. The experimental values were presented in Ref. 5. The units are in kcal/mol.

ATP

Pyrophosphate

Hydrolysis

ΔGhyd (PW)

ΔGhyd (Expt.)

R1

−6.6



R2

−25.2

−9.9

R3

−16.3

−10.7

R4

−14.3

−9.5

R5

−7.3

−7.5

R6

−17.8

−7.7

R7

−16.2

−7.1

R8

−17.8

−10.4

20

ACS Paragon Plus Environment

5

Page 21 of 29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure captions: Figure 1

Snapshot of the QM/MM-ER simulation for ATP4− embedded in water droplet.

The blue surface represents the iso-surface at 0.02 a.u. of electron density n! of ATP4−. Figures 2

(a) Energy distribution functions ρ(ε) in solution and ρ0(ε) in pure solvent

systems for the calculation of the free energy Gsol(ATP4−). In the inset, snapshots of the water molecules taken at the intervals of 50 fs during a 550-fs simulation are superimposed. The water molecules with red oxygen atoms have the pair interactions ε with the solute in the region of −35 < ε < −30 kcal/mol, while those with blue ones lie in −30 < ε < −25 kcal/mol. (b) Energy distribution functions ρ(ε) in solution and ρ0(ε) in pure solvent for the calculation of Gsol(ADP3−). Figures 3

Energy distribution functions P(η) in solution and P0(η) in reference systems for

the calculations of the electronic free energies (a) Gele(ATP4−) and (b) Gele(ADP3−). The values of η are actually shown relative to the average in the solution system of interest. Figure 4

Free energy changes ΔG (ΔGhyd, ΔGsol, and ΔGele) for the hydrolyses of ATP and

PPi with respect to the variation of numbers Ne of the excess electrons. The explicit reaction schemes are presented in Table 1. The free-energy values at Ne = 0 for the ATP and PPi systems are so close that the points for ATP are hardly discernible within the resolution of the figure. Figure 5

Solvent effect in the free energy ΔGsol for ATP hydrolyses with respect to the

variation of numbers Ne of the excess electrons for QM/MM calculations in solvent water and classical MD calculations in solvent water and ethanol.

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Contents (TOC)

22

ACS Paragon Plus Environment

Page 22 of 29

Page 23 of 29

Figure 1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ACS Paragon Plus Environment

The Journal of Physical Chemistry

Distribution / mol/kcal

−30

ATP 4−

−35

Solution

−20

Pure

−25 −15 −10 −5 0 5 10 ACS Paragon Plus Environment

40

35

30

25

20

15

10

5

0 −40 Energy coordinate / kcal/mol

Figure 2(a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 29

Page 25 of 29

Distribution / mol/kcal

−30

ADP 3−

−35

Solution

−20

Pure

−25 −15 −10 −5 0 5 10 ACS Paragon Plus Environment

40

35

30

25

20

15

10

5

0 −40 Energy coordinate / kcal/mol

Figure 2(b)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The Journal of Physical Chemistry

Distribution / mol/kcal

ATP 4−

−25

Solution

Reference

−20 −15 −10 −5 ACS Paragon Plus Environment

0.40

0.35

0.30

0.25

0.20

0.15

0.10

0.05

0 −30 Energy coordinate / kcal/mol

Figure 3(a)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 29

Page 27 of 29

Distribution / mol/kcal

ADP 3−

−25

Solution

Reference

−20 −15 −10 −5 ACS Paragon Plus Environment

0.40

0.35

0.30

0.25

0.20

0.15

0.10

0.05

0 −30 Energy coordinate / kcal/mol

Figure 3(b)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

The Journal of Physical Chemistry

∆G (kcal/mol) 300

200

100

0

−100

−200

−300

−400 0

∆Ghyd (ATP)

∆Gsol (ATP)

3 4 ACS Paragon Plus Environment

∆Gele (ATP)

∆Ghyd (PPi)

∆Gsol (PPi)

∆Gele (PPi)

2

Born (ε = 78.0)

1

Number of excess electrons Ne

Figure 4

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 29

Page 29 of 29

∆Gsol (kcal/mol) 150

100

50

0

QM/MM in water classical MD in water classical MD in ethanol

0 3 4 Number of excess electrons Ne

Figure 5

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

ACS Paragon Plus Environment