Ab Initio Molecular Dynamics Simulations of Amino ... - ACS Publications

Subscriber access provided by UNIV OF CAMBRIDGE

Article

Ab initio Molecular Dynamics Simulations of Amino Acids in Aqueous Solutions: Estimating pKa Values from Metadynamics Sampling Anil Kumar Tummanapelli, and Sukumaran Vasudevan J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.5b05211 • Publication Date (Web): 25 Aug 2015 Downloaded from http://pubs.acs.org on August 29, 2015

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 26

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

Ab initio Molecular Dynamics Simulations of Amino Acids in Aqueous Solutions: Estimating pKa Values from Metadynamics Sampling

Anil Kumar Tummanapelli and Sukumaran Vasudevan* Department of Inorganic and Physical Chemistry Indian Institute of Science, Bangalore 560012, INDIA

* Author to whom correspondence may be addressed. E-mail: [email protected]. Tel: +91-80-2293-2661. Fax: +91-80-2360-1552/0683;

ACS Paragon Plus Environment

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

Page 2 of 26

ABSTRACT Changes in protonation-deprotonation of amino acid residues in proteins play a key role in many biological processes and pathways. Here we report calculations of the free energy profile for the protonation-deprotonation reaction of the 20 canonical α-amino acids in aqueous solutions using ab initio Car-Parrinello molecular dynamics simulations coupled with metadynamics sampling. We show here that the calculated change in free energy of the dissociation reaction provides estimates of the multiple pKa values of the amino acids that are in good agreement with experiment. We use the bond length dependent number of the protons coordinated to the hydroxyl oxygen of the carboxylic and the amine groups as the collective variables to explore the free energy profiles of the Bronsted acid-base chemistry of amino acids in aqueous solutions. We ensure that the amino acid undergoing dissociation is solvated by at least three hydrations shells with all water molecules included in the simulations. The method works equally well for amino acids with neutral, acidic and basic side chains and provides estimates of the multiple pKa values with a mean relative error, with respect to experimental results, of 0.2 pKa units.


2

Page 3 of 26

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


INTRODUCTION Proton transfer plays a critical role in the function of proteins and peptides and is one of the most frequently occurring reactions in nature. The pH-dependent changes in protonationdeprotonation of amino acid residues in proteins play a key role in many biological processes and pathways–from enzyme activity to energy conversion in living cells.

Quantitative

characterization of such processes requires knowledge of the acid equilibrium constant, Ka, of the molecules involved; it is this quantity that determines the state of protonation of the system at a given pH. It is, therefore not surprising that the theoretical estimation of the pKa (= - log Ka) values of the titratable groups in proteins has been considered as one of the most important and challenging problems in biophysical chemistry.1 As the first step in this direction is the need for procedures, free of empirical scaling or corrections, that can provide accurate estimates of the values of the pKa s’ of the 20 canonical α- amino acids. Predicting pKa values for acid/base reactions in aqueous solutions is one of the key challenges in theoretical chemistry. This has been an active area of research with focus on developing new methods and procedures that can provide estimates of pKa values that are comparable, within 0.5 pKa units, with experiment and applicable for a variety of different chemical systems.1-3 This is a demanding exercise for an accuracy of 0.74 pKa units requires an errors of less than 1 kcal/mole in estimates of the free energy for the dissociation reaction.4 Most methods for estimating acid dissociation constants involve constructing a suitable thermo-chemical cycle with the dissociation free energy evaluated from a summation of the free energy change associated with the individual steps. These steps include the gas-phase dissociation free energy


3


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 26

for the acid molecule and the solvation free energy difference of the neutral acid and the dissociated products. Changes in free-energy for the gas-phase reaction can be evaluated with sufficient accuracy but in solution, where hydration effects have to be considered, the same is not true. The most difficult part in the theoretical prediction of pKa values is to correctly account for the structure and dynamics of the solvent molecules solvating the acid molecule. This is both conceptually and computationally demanding. In most calculations the dissociation of the solute molecule is handled with rigor, at the quantum mechanical (QM) level , but the solvent modeled as a structure less dielectric continuum.5 A major shortcoming of the continuum solvent model is that short-range intermolecular interactions, that are particularly important while modeling hydrated ionic species are not explicitly incorporated in the calculations. The short range interactions include hydrogen bonding and ion-dipole interactions between solvent molecules. Specific inclusion of a finite number of water molecules in the dissociation step is one of the approaches to enhance the accuracy of the calculations.6-8 Explicit solvation shell models have been used to estimate the pK’s of select amino acids, in but agreement with experimental values is not too satisfactory.9 An alternate approach is to model the solvent molecules by force-field based molecular mechanics methods (MM) while describing the dissociation process by quantum chemical techniques.10-13 It has been shown that boundary issues inherent to the QM-MM method can be avoided using reactive-MD procedure that are appropriately parameterized. Multiscale reactive–MD simulations have been used to predict absolute pKa values of aspartic and glutamic amino acids.14 A critical factor in any dissociation event is to account for the structure and dynamics of water molecules hydrating the neutral and charged dissociating species. This is best described by ab initio or quantum MD simulations that consider all solvent molecules explicitly at the same level


4

Page 5 of 26

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 theory.15 In these simulations density functional theory (DFT) is used to describe the electronic structure of both the dissociating acid and the solvating water molecules. Short range interactions involving solvent molecules are explicitly accounted for in the simulations, unlike in continuum models. Increased access to computer resources have made these methods popular and a number of investigations of acid dissociation reactions using ab initio MD simulations including Car-Parrinello molecular dynamic(CPMD), have been reported.16-19 The acid molecules investigated were HCOOH, CH3COOH and H2CO3.16-18 These simulations do not directly provide pKa values but procedures have been developed to obtain the acid dissociations constants. The dissociation of histidine in aqueous media

had

been investigated using

constrained CPMD with free energies obtained by integrating the potentials of mean force with dissociation constants obtained after subtracting profile for the self-dissociation of water.19 Dissociation free energies have also been calculated using the distribution of vertical energy gaps for the addition or removal of a proton obtained by combining DFT based MD simulations with thermodynamic integration.20 The procedure has been used with a slight modification to obtain acidity constants for select set of amino acids. The computed pK values showed a mean error of 2.1 pK units with respect to experimental values.21 The dissociation of a weak acid is a rare event and consequently exceptionally long times are required before one is observed in normal simulation procedures. Metadynamics sampling is one of the successful techniques that can handle rare and infrequent events in MD simulations.22-24 The formalism prevents the system from being trapped in a free-energy minima, allowing it to explore the entire free energy space without revisiting regions where it has been in the past .It does so by biasing the dynamics with a history dependent potential, modeled by repulsive inverted Gaussians, which

are subsequently

dropped during propagation. The bias acts on


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

Page 6 of 26

select degrees of freedom, referred to as collective variables. In the present study these collective variables, depending on the reaction, are either the distance-dependent number of protons coordinating the carboxylic or

amine groups of the amino acid. Metadynamics sampling is

relatively quick and efficient and , therefore, ideally suited for studying dissociation of weak acids.15-18,25,26 The method has been used successfully to investigate the dissociation of acetic acid in water.17,25 Minima associated with the neutral and ioinized states of the acid along with a shallow minimum corresponding to contact ion pair formation could clearly be distinguished in the metadynamics computed free energy profile. The pKa value estimated from the difference in the values of the free energy minima corresponding to the neutral and dissociated acetic acid was in good agreement with experiment.25 The formalism was shown to correctly account for the influence of the inductive effect as well as hydrogen bonding on pKa values of weak organic acids in aqueous solutions and also provide accurate estimates of the successive dissociation constants of polyprotic acids in water.25,26 A notable highlight of the

combined CPMD-

metadynamics sampling procedure is that it retains conceptual simplicity while still considering all solvent water molecules explicitly in the simulation. As compared to simple organic weak acids the amino acids pose the additional challenge that they are zwitterions by the fact that the –NH2 group is a base while the –COOH an acid. In addition, they may possess side-groups, which could be non-polar or have a dipole moment or be either acidic or basic in nature. The task therefore is to be able to estimate the multiple pKa s that characterizes the acid-base chemistry of the amino acids in aqueous solutions. Here we show that CPMD27 simulations in conjunction with metadynamics computation of the free-energy landscape of the protonation–deprotonation reaction can provide accurate estimates of the multiple pKa’s of the 20 canonical amino acids in aqueous solution.


6

Page 7 of 26

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


SIMULATION METHODOLOGY Ab initio NVT MD simulations and the free energy computation using the extended Lagrangian version of metadynamics were performed using the CPMD v3.11.2 program28 Simulations were initiated from an amino acid molecule with water molecules, sufficient to complete three hydration shells, placed in a cubic box with periodic boundary conditions. A density of unity for the solution was achieved by adjusting the volume of the simulation cell so that the space available for the water molecules, after subtracting the van der Waals volume of the amino acid, corresponds to a water density of unity. The number of water molecules and the dimension of the simulation cell are different for each amino acid and are specified in the Supporting Information, Table S1. The temperature of the simulation cell was maintained at 300K by a Nose-Hoover thermostat. The DFT-HCTH/120 functional with the empirical dispersion correction of Grimme was used.29,31 Wave functions were optimized till a convergence value of 5 × 10-8 a.u.. Details of the simulation methodology are similar to that reported earlier.25,26 In the metadynamics computation of the free-energy for dissociation process the collective variable used was the bond-distance dependent coordination number of protons bound to the dissociating group , nXH (X= O and N), defined as

∑

(

)

(

)

where is rXH is the X-H bond distance and rc (= 1.6 Å) the critical distance beyond which the X-H bond breaks The summation runs over protons of all water molecules in the cell and that of the dissociating group. The exponents ensures a smooth decay of

nXH.23 For protonation-

deprotonation of the carboxylic group the collective variable nOH has a value close to unity in the


7


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 8 of 26

un-dissociated state and is zero on dissociation. For the amine group the collective variable nNH defines the state of protonation of the –NH2 group and has a value close to unity in the protonated –NH3+ state and is zero on deprotonation. In the present study the simulations were initiated from the amino acid in the fully protonated state i.e. both nOH and nNH have a value of unity, and successively deprotonated. This, however, is a matter of convenience, simulations initiated with the same amino acid in different degrees of protonation gave similar results. In these calculations the simulation cell remained charge neutral at all instances; thus when the amino acid was fully protonated electrical neutrality was preserved by the presence of a negatively charged hydroxyl ion in the cell. After the dissociation of the carboxylic group the product hydronium ion remains in the cell to ensure charge neutrality during simulations of the deprotonation of the amine group. The free energy profile was sampled by biasing the dynamical time evolution path with Gaussian bias potentials of height 0.0005 a.u. and width 0.1 a.u. that were deposited at time intervals of 0.04 ps. The Gaussian height was chosen after ascertaining that calculations using a smaller Gaussian height, 0.0001 a.u., gave similar results. The sampling was continued till the entire free-energy landscape was traversed and the motion in the collective variable space was diffusive. For simulations of the dissociation of amino acids with acidic side groups nOH is defined for a designated carboxylic group. These simulations, too, were initiated from the fully protonated amino acid and successively de-protonated. The hydronium ion formed on dissociation of the first carboxylic group remains as part of the simulation cell to ensure charge neutrality for the computation of the dissociation constant of the second carboxylic group. Similarly for the protonation-deprotonation of the amine groups of the amino acids with basic side-groups the


8

Page 9 of 26

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


collective variable nNH is defined for a designated amine group. Here, too, hydroxyl anions were present in the cell to ensure charge neutrality during simulations of the deprotonation of the second amine group. It may be emphasized that in these calculations the results were the same irrespective of whether the simulations were started from the fully protonated amino acid and was sequentially de-protonated, or if initiated from the de-protonated state and successively protonated. RESULTS AND DISCUSSIONS The protonation–deprotonation reactions and the associated acid equilibrium constants of amino acids with neutral side-chains, acidic side-chains and basic side-chains are discussed separately in the following sections. Amino acids with Neutral Side-chains. The simulations were initiated from the fully protonated positively charged state (Scheme-1) and successively deprotonated. The definitions of the pKas are shown in Scheme-1. In the following discussion we consider the protonation–deprotonation of glycine (R = H) as a prototype for this class of amino acids. The free energy profile for the pKa1 acid dissociation process as defined by the collective variable nOH, is shown in Figure 1a. These calculations involved placing a protonated glycine in a cubic box of side 12.0Å along with 55 water molecules and an OH- hydroxyl anion. The number of water molecules ensured that at least three hydration shells were complete and the dimensions of the box that the water density as

Scheme-1


9


0 2 3

-2

Free Energy, kcal/mol

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 10 of 26

∆G = 3.31 kcal/mol pKa1 = 2.41(Exptl. 2.35)

-4 1 -6 0

2'

3'

-4

∆G = 13.54 kcal/mol pKa2 = 9.86 (Exptl. 9.78)

-8

-12 -16 -20 1.0

1' 0.7

0.4

0.1

Collective Variable, nXH

Figure 1 (a) The free-energy profiles for the dissociation of the carboxylic group of glycine. The difference in free energies of the neutral and ionized species and the estimated pKa1 value is indicated (experimental values are in parentheses). (b) Snapshots along the trajectory for the carboxylic group dissociation of glycine at the locations indicated in Figure 1a. The proton being abstracted is shown in yellow and the water molecules involved in the dissociation reaction highlighted. (c) The free-energy profiles for the deprotonation of the NH3+, group of glycine. The difference in free energies of the protonated and deprotonated species, as well as the estimated pKa2 value is indicated(experimental values are in parentheses). (d) Snapshots along the trajectory for the deprotonation of the NH3+, group of glycine. The snapshots are for the locations indicated in Figure 1c. The proton being abstracted is shown in yellow and the water molecules involved in the dissociation reaction highlighted. Movie files of configurations along the trajectory for the dissociation of the carboxylic and the protonation of the amine are provided in the Supporting Information, Movie S1 and S2).


10

Page 11 of 26

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


close to unity. During the metadynamics simulations complete dissociation of the carboxylic group is observed within 6 ps (see Supporting Information, Figure S1a). At this stage the system escapes from the product (dissociated) state and the motion is diffusive and multiple dissociation (nOH ≈ 0) and association (nOH ≈ 1) events are observed (see Supporting Information, Figure S1a). In Figure 1a the minima at nOH ≈1 corresponds to the neutral un-dissociated carboxylic group while that at nOH ≈ 0 to the dissociated ionized –COO-. Snapshots (Figures 1b) at the locations indicated on Figure 1a confirm that these minima correspond to the neutral and ionized species, respectively. 2 ps after the metadynamics simulations are initiated the carboxylic proton of glycine is abstracted by a water molecule but the hydronium ion and the resultant carboxylate ion remain as a bound pair that manifests as a shallow minima, labeled 2, in the free energy profile (Figure 1a). The dissociation process is completed with the break-up of the contact-ion pair . The free-energy profile in Figure 1a are identical irrespective of whether the simulations were initiated from the neutral carboxylic or dissociated carboxylate species (see Supporting Information, Figure S2). Free energy profiles computed with a larger simulation cell containing 100 water molecules were identical to that shown in Figure 1a (see Supporting Information, Figure S3) and consequently subsequent simulations were done with 55 water molecules in the simulation cell. Subsequent to the dissociation of the carboxylic group the deprotonation of the NH3+ group of glycine was investigated. The free energy profile for the pKa2 deprotonation process, as defined by the collective variable nNH, is shown in Figure 1b. In these simulations the deprotonation of the NH3+ is observed within the first 6.2 ps (see Figure S1b). The free energy profile shows two well defined minima corresponding to the protonated and deprotonated states of the amine group of glycine. The free energy profiles for simulations initiated from either the protonated or de-


11


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 26

protonated states of the glycine amine group were identical (see Supporting Information, Figure S4). The free energy landscape computed by metadynamics simulations of aqueous solutions of glycine using the distance-dependent coordination numbers, nOH and nNH, as collective variables (Figure 1a and Figure 1c) are clearly able to distinguish the neutral and ionized states of the carboxylic and amine groups of glycine as individual well-defined minima . The dissociation constants may be estimated from the free energy profiles in Figure 1a and 1c as the difference in the values of the free energy at the minima corresponding to the protonated and deprotonated species, pKa = ∆G / 2.303 RT. The estimated value of pKa1, the carboxylic dissociation is 2.41 while the value of pKa2, the deprotonation of the amine NH3 is 9.86. These values are in good agreement with the experimental values of pKa1 and pKa2 for glycine, 2.35 and 9.78, respectively. The estimated value of the isoelectric point, pI, the pH where the glycine molecule exists in an uncharged state, obtained as the average of pKa1 and pKa2, is 6.14 ; the experimentally determined value, 5.97. The estimated and experimental pKa values of glycine are in good agreement in spite of the many approximations used in this computation. These approximations include the choice of the functional, ignoring quantum effects as well as zero-point energies32 and, in addition, there are statistical errors associated with the metadynamics sampling technique. It had, however, been pointed out in earlier reports that since the pKa values are estimated from the difference in free energy values associated with the protonated and deprotonated species, there is a cancellation of errors.25,26 It is for this reason that the procedure is able to accurately predict the pKa values of polyprotic organic acids

and correctly account for the influence of the inductive effect and

hydrogen bonding on acid dissociation constants.25,26 It is, therefore, not surprising that the


12

Page 13 of 26

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


estimate of the pKa values of glycine, obtained as the difference in free-energy values leading to a cancellation of errors, are in good agreement with experiment. In the following sections we use the CPMD-metadynamics formalism to estimate the pKa values of the remaining canonical amino acids in aqueous solution. The protonation–deprotonation of amino acids with neutral side-chains was investigated by procedures similar to that for glycine. The number of water molecules and the dimensions of the simulation cell were different for each of the amino acids (see Supporting Information, Table S1). The number of water molecules and the dimension of the cell were chosen so as to ensure that at least three hydration shells around the amino acid were complete and that the density of water as close to unity. The free energy profiles of the dissociation of the carboxylic group, as defined by the collective variable nOH, were computed by the metadynamics simulations and so also the deprotonation of the amine, NH3+, using nNH as the collective variable. For each of the amino acids the free-energy profiles showed two distinct minima corresponding to the different states of protonation (Supporting Information, Figures S5-S18). The pKa1 and pKa2 values estimated from the difference in the values of the free-energy minima, corresponding to the protonated and deprotonated states, are tabulated in Table-1 along with the experimentally determined values of these quantities. The estimated and experimental values of the isoelectric point, pI, obtained from the average of pKa1 and pKa2 , are also included in Table-1. The agreement between the estimated and experimental values of the pKas is good. Amino acids with Acid Side-chains. The simulations were initiated from the fully protonated state of the amino acid (Scheme-2) and successively deprotonated. It may be emphasized that the


13


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 14 of 26

Scheme-2

Figure 2: The free energy profiles for the deprotonation of aspartic acid. For the curves, I and II, the deprotonation of the –COOH

groups the collective variable is nOH while for the

deprotonation of the NH3+ groups (curve III) the collective variable is nNH. choice of the initial state of protonation is a matter of convenience; the results are identical if instead of proceeding from right to left in Scheme 2, the protonation process, proceeding from left to right in Scheme 2, was investigated (see Supporting Information, Figure S2). The definition of the pKa s of amino acids with acid side chains is indicated in Scheme-2. The freeenergy profiles for the three deprotonation reactions of the fully protonated aspartic acid are


14

Page 15 of 26

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


shown in Figure 2. The profiles are for the dissociation of the amino acid carboxylic group followed by the dissociation of the side chain carboxylic group and finally the deprotonation of the NH3+, group. For the first two processes the collective variable is nOH, the distance-dependent coordination number of protons attached to the carboxylic group undergoing dissociation while for the deprotonation of the NH3+ it is nNH. For each of the free energy profiles two distinct minima corresponding to the protonated and deprotonated states are observed and the value of the pKa s estimated from the difference in the free energy values at the minima (Figure 2). The shallow minima (Figure 2) in the free energy profiles correspond to a contact-ion pair that subsequently breaks-up to complete the dissociation process. For aspartic acid the value of pKa1, pKaR and pKa2 , 2.06, 4.0 and 9.99, respectively are in good agreement with the experimentally determined values, 1.99, 3.9 and 9.9, respectively. For glutamic acid, too, the free energy profiles for each of the three deprotonation reactions showed two distinct minima corresponding to the protonated and deprotonated states (Supporting Information, Figure S19) and the estimates of the pKa values in good agreement with experiment (Table-1). Amino acids with Basic Side-chains. The definition of the pKas for amino acids with basic side chains is indicated in Scheme-3. For lysine and arginine pKaR > pKa2 and the sequence of

Scheme-3


15


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 16 of 26

Figure 3: The free energy profiles for the deprotonation of histidine. For the curve I, the deprotonation of the -COOH the collective variable is nOH while for the deprotonation of the NH3+ groups (curves II and III) the collective variable is nNH. deprotonation steps in the simulations are reversed from that shown in Scheme-3. The free energy profiles for the dissociation of the carboxylic group and the NH3+ groups of histidine are shown in Figure 3 (the free energy profiles for the deprotonation steps of lysine and arginine are provided as part of the Supporting Information, Figures S20,S21). The definition of the collective variables nXH (X =O,N) are as discusssed in the previous sections. Each of the profiles show two well-defined minima, corresponding to the different states of protonation, and the pKa values estimated from the difference in the values of the free energy at


16

Page 17 of 26

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 1: Values of pKa1, pKa2, pKaR and pI estimated from the metadynamic simulations for the 20 canonical amino acids in aqueous solutions. The values in parenthesis are the experimentally determined values.33,34

Side Group

Neutral

Amino Acids

pKa1

pKa2

pKaR

pI

Glycine, Gly

2.41(2.35)

9.86(9.78)

-

6.14(5.97)

Alanine, Ala

2.42(2.35)

9.96(9.87)

-

6.19(6.02)

Valine, Val

2.37(2.29)

9.84(9.74)

-

6.11(5.97)

Leucine, Leu

2.41(2.33)

9.85(9.74)

-

6.13(5.98)

Isoleucine, Ile

2.39(2.32)

9.86(9.76)

-

6.13(6.02)

Phenylalanine, Phe

2.26(2.20)

9.40(9.31)

--

5.83(5.48)

Tryptophan, Trp

2.52(2.46)

9.49(9.41)

-

6.01(5.88)

Serine, Ser

2.27(2.19)

9.33(9.21)

-

5.80(5.68)

Threonine, Thr

2.16(2.09)

9.22(9.10)

-

5.69(6.53)

Methionine, Met

2.19(2.13)

9.39(9.28)

-

5.79(5.75)

Asparagine, Asn

2.21(2.14)

8.91(8.72)

-

5.56(5.41)

Glutamine, Gln

2.25(2.17)

9.23(9.13)

-

5.74(5.65)

Proline, Pro

2.13(1.95)

10.74(10.64)

-

6.44(6.10)

Tyrosine, Tyr

2.27(2.20)

9.30(9.21)

-

5.79(5.65)

Cysteine, Cys

2.00(1.92)

10.80(10.70)

-

5.24(5.14)

Aspartic Acid, Asp

2.06(1.99)

9.99(9.90)

4.00(3.90)

3.03(2.87)

Glutamic Acid, Glu

2.16(2.10)

9.55(9.47)

4.15(4.07)

3.16(3.22)

Histidine, His

1.89(1.80)

9.43(9.33)

6.13(6.04)

7.78(7.58)

Lysine, Lys

2.23(2.16)

9.17(9.06)

10.63(10.54)

9.90(9.74)

Arginine, Arg

1.90(1.82)

9.11(8.99)

12.58(12.48)

10.85(10.76)

Acidic

Basic


17


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 18 of 26

Figure 4: Graphical comparison of the estimated and experimentally determined pKa1, pKa2 and pKaR values of the 20 amino acids (experimental values are from Ref. 33 and 34). the minima. The estimated and experimentally determined pKa values for the amino acids with basic side chains are in good agreement with corresponding experimentally determined values (Table 1). The results of the pKa values of amino acids with neutral, acid and basic side chains estimated using the CPMD simulations in conjunction with metadynamics computations of the free energy profiles for the protonation-deprotonation reaction are summarized in Table-1 and Figure 4. For comparison the experimentally determined values of these quantities are also included in the table and figure. It may be seen that the estimated and experimental values of the pKas of the


18

Page 19 of 26

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


amino acids are in very good agreement. The mean error with respect to experiemntal results is 0.2 pKa units with a maximum error of 0.48 pKa units. It is important to compare the results of this study on estimating

pKa values of amino acids using CPMD in conjuction with

metadynamics simulation of the free energy with previous ab initio studies. The pKa values of a select number of amino acids had been reported by constructing the appropriate thermochemical cycle with solvation free energies computed using the explicit solvation shell model. The mean error with respect to experimental results for this study is 3.0 pKa units.9 As mentioned in the introduction the estimates of pK values by ab initio molecular dynamics along with thermodynamic integration showed a mean error of 2.1 pK units with respect to experimental values.21 The results of the present study show a much lower value of the mean error as compared to earlier reports.9,21 It should also be noted that in the present study the mean error of 0.2 pKa units was computed considering the multiple values of the pKa s for all 20 canonical amino acids. CONCLUSIONS In conclusion we have shown that CPMD based metadynamics simulations can predict pKa values of amino acids with good accuracy. The procedure outlined here considers solvent water molecules explicitly, is free of any empirical corrections or scaling parameters and is conceptually simple. For each of the amino acids the number of water molecules considered was sufficient to complete the three hydration shells. We use the bond-length dependent number of the protons bound to the hydroxyl oxygen of the carboxylic group as the collective variable to explore the free energy profile of the dissociation of the carboxylic groups and the bond-length dependent number of protons attached to the amine as the collective variable to define the free energy profile of the protonation-deprotonation of the amine group. Two distinct minima


19


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 20 of 26

corresponding to the protonated and deprotonated states for each of the acid groups of the amino acid are observed and the pKa values estimated from the difference in their free energy values at the minima. Since the pKa values are estimated as the difference in free-energy values the procedure benefits from a reduction in errors due to cancellation, which appears to hold true for all twenty of the canonical amino acids. The agreement between the ab initio estimated and experimentally determined values of the multiple pKa s

for all 20 canonical amino acids investigated is good. The values of the mean

error with respect to experiemntal results, 0.2 pKa units, and the maximum error, 0.48 pKa units, is well within the desirable prediction accuracy for pKa values. We believe that the present study lays the groundwork for future ab initio studies that would be able to provide accurate estimates for the acid equilibrium constants of simple peptides, and eventually to predict the pKa values of amino acid residues in proteins.

ASSOCIATED CONTENT Supporting Information Available: (Table S1) The number of water molecules and dimension of the simulation cell for the amino acid molecules investigated. (Figure S1) Evolution of the collective variables, nOH and nNH ,during a 10ps metadynamics simulation of the dissociation of the carboxylic and amine groups of glycine. (Movie S1) Movie file of the configurations along the trajectory of the deprotonation of the carboxylic group of glycine. (Movie S2) Movie file of the configurations along the trajectory of the deprotonation of the NH3+ group of glycine. (Figure S2) The free–energy profiles and trajectories for the dissociation and association of the carboxylic group of glycine. (Figure S3) The free-energy profiles for the dissociation of the carboxylic group of glycine from simulations with 55 water molecules and 100 water molecules.


20

Page 21 of 26

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


(Figure S4) The free–energy profiles and trajectories for the deprotonation and protonation of the amine group of glycine. (Figures S5-S21) The free energy profiles for the protonation– deprotonation of a series of amino acids in aqueous solutions. This material is available free of charge via the Internet at http://pubs.acs.org.

ACKNOWLEDGMENTS The authors thank the Supercomputer Education and Research Centre at the Indian Institute of Science, Bangalore for use of the IBM Blue Gene cluster computational facility.

REFERENCES (1) Alexov, E.; Mehler, E. L.; Baker, N.; Baptista, A. M.; Huang, Y.; Milletti, F.; Nielsen, J. E.; Farrell, D.; Carstensen, T.; Olsson, M. H.; Shen, J. K.; Warwicker, J.; Williams, S.; Word, J. M. Progress in the Prediction of pKa Values in Proteins. Proteins 2011, 79, 3260-3275. (2) Alongi, K. S.; Shields, G. C. Theoretical Calculations of Acid Dissociation Constants: A Review Article. Annu. Rep. Comput. Chem. 2010, 6, 113-138. (3) Ho, J.; Coote, M. L. A Universal Approach for Continuum Solvent pKa Calculations: Are we there Yet? Theor. Chem. Acc. 2010, 125, 3-21. (4) Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. Absolute pKa Determinations for Substituted Phenols. J. Am. Chem. Soc. 2002, 124, 6421-6427.


21


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 22 of 26

(5) Lim, C.; Bashford, D.; Karplus, M. Absolute pKa Calculations with Continuum Dielectric Methods. J. Phys. Chem. 1991, 95, 5610-5620. (6) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Adding Explicit Solvent Molecules to Continuum Solvent Calculations for the Calculation of Aqueous Acid Dissociation Constants. J. Phys. Chem. A 2006, 110, 2493-2499. (7) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Aqueous Solvation Free Energies of Ions and Ion−Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton. J. Phys. Chem. B 2006, 110, 16066-16081. (8) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. SM6: A Density Functional Theory Continuum Solvation Model for Calculating Aqueous Solvation Free Energies of Neutrals, Ions, and Solute−Water Clusters. J. Chem. Theory Comput. 2005, 1, 11331152. (9) Gupta, M.; Silva, E. F. d.; Svendsen H. F. Explicit Solvation Shell Model and Continuum Solvation Models for Solvation Energy and pKa Determination of Amino Acids. J. Chem. Theory Comput. 2013, 9, 5021-5037. (10) Choi, C. H.; Re, S.; Feig, M.; Sugita, Y. Quantum Mechanical/Effective Fragment Potential Molecular Dynamics (QM/EFP-MD) Study on Intra-Molecular Proton Transfer of Glycine in Water. Chem. Phys. Lett. 2012, 539, 218−221. (11) Ghosh, M. K.; Re, S.; Feig, M.; Sugita, Y.; Choi, C. H. Interionic Hydration Structures of NaCl in Aqueous Solution: A Combined Study of Quantum Mechanical Cluster Calculations and QM/EFPMD Simulations. J. Phys. Chem. B 2013, 117, 289−295.


22

Page 23 of 26

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


(12) Ghosh, M. K.; Uddin, N.; Choi, C. H. Hydrophobic and Hydrophilic Associations of a Methanol Pair in Aqueous Solution. J. Phys. Chem. B 2012, 116, 14254−14260. (13) Uddin, N.; Choi, T. H.; Choi, C. H. Direct Absolute pKa Predictions and Proton Transfer Mechanisms of Small Molecules in Aqueous Solution by QM/MM-MD. J. Phys. Chem. B 2013, 117, 6269−6275. (14) Nelson, J. G.; Peng, Y.; Silverstein, D. W.; Swanson, J. M. J. Multiscale Reactive Molecular Dynamics for Absolute pKa Predictions and Amino Acid deprotonation. J. Chem. Theory Comput. 2014, 10, 2729−2737. (15) Hassanali, A. A.; Cuny, J.; Verdolino, V.; Parrinello, M. Aqueous Solutions: State of the Art in Ab initio Molecular Dynamics. Philos. Trans. R. Soc. A 2014, 372, 20120482. (16) Lee, J. G.; Asciutto, E.; Babin, V.; Sagui, C.; Darden, T.; Roland, C. Deprotonation of Solvated Formic Acid: Car-Parrinello and Metadynamics Simulations. J. Phys. Chem. B. 2006, 110, 2325-2331. (17) Park, J.M.; Laio, A.; Iannuzzi, M.; Parrinello, M. Dissociation Mechanism of Acetic Acid in Water. J. Am. Chem. Soc. 2006, 128, 11318-11319. (18) Galib, M.; Hanna, G. Mechanistic Insights into the Dissociation and Decomposition of Carbonic Acid in Water via the Hydroxide Route: An Ab Initio Metadynamics Study. J. Phys. Chem. B. 2011, 115, 15024-15035. (19) Ivanov, I.; Chen, B.; Raugei, S.; Klein, M. L. Relative pKa Values from First-Principles Molecular Dynamics: The Case of Histidine Deprotonation. J. Phys. Chem. B 2006, 110, 6365−6371.


23


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 26

(20) Sulpizi, M.; Sprik, M. Acidity Constants from Vertical Energy Gaps: Density Functional Theory Based Molecular Dynamics Implementation. Phys. Chem. Chem. Phys. 2008, 10, 5238–5249. (21) Mangold, M.; Rolland, L.; Costanzo, F.; Sprik, M.; Sulpizi, M.; Blumberger, J. Absolute pKa Values and Solvation Structure of Amino Acids from Density Functional Based Molecular Dynamics Simulation. J. Chem. Theory Comput. 2011, 7, 1951. (22) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. PNAS 2002, 99, 12562-12566. (23) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2003, 90, 238302-4. (24) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 826. (25) Tummanapelli, A. K.; Vasudevan, S. Dissociation Constants of Weak Acids from Ab Initio Molecular Dynamics Using Metadynamics: Influence of the Inductive Effect and Hydrogen Bonding on pKa Values, J. Phys. Chem. B, 2014, 118, 13651–13657. (26) Tummanapelli, A. K.; Vasudevan, S. Estimating successive pKa values of polyprotic acids from ab initio molecular dynamics using metadynamics: the dissociation of phthalic acid and its isomers. Phys. Chem. Chem. Phys., 2015, 17, 6383-6688. (27) Car, R.; Parrinello, M.

Unified Approach for Molecular Dynamics and Density-

Functional Theory. Phys. Rev. Lett. 1985, 55, 2471-2474. (28) CPMD, version 3.13.2; IBM Corp 1990−2008, MPI für Festkörper-forschung Stuttgart 1997−2001. http://www.cpmd.org/.

(29) Boese, A. D.; Doltsinis, N. L.; Handy, N. C.; Sprik, M. New Generalized Gradient Approximation Functionals. J. Chem. Phys. 2000, 112, 1670-1678.


24

Page 25 of 26

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


(30) S. Izvekov and G. A. Voth, Ab Initio Molecular Dynamics Simulation of Aqueous Proton Solvation and Transport Revisited, J. Chem. Phys., 2005, 123, 044505-9. (31) S. Grimme, Semiempirical GGA-Type Density Functional Constructed with a LongRange Dispersion Correction, J. Comput. Chem., 2006, 27, 1787−1799. (32) M. Garcia-Viloca, C. Alhambra, D. G. Truhlar and J. Gao, Inclusion of QuantumMechanical Vibrational Energy in Reactive Potentials of Mean Force, J. Chem. Phys., 2001, 114, 9953-9958. (33) Dawson, R. M. C.; Elliott, D. C.; Elliott, W. H.; Jones, K. M. Data for Biochemical Research, 3rd ed.; Clarendon Press: Oxford, 1986. (34) Edsall, J. T.; Wyman, J. Biophysical Chemistry; Academic Press: New York, 1958.


25


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 26

Table of Content

pKa1 2.41 2.35

pKa1 +

NH3-CH2-COOH

pKa2 +

-

NH3-CH2-COO

-

NH2-CH2-COO

9.86

pKa2 9.78


26

Ab Initio Molecular Dynamics Simulations of Amino ... - ACS Publications

Recommend Documents