Correction to Constant-pH Molecular Dynamics Simulations for Large

If it takes too long to load the home page, tap on the button below. Try Again! .... Cancel. To get a pairing code: ... Find the code on the page and ...
0 downloads 0 Views 243KB Size
Erratum pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Correction to Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems Brian K. Radak, Christophe Chipot, Donghyuk Suh, Sunhwan Jo, Wei Jiang, James C. Phillips, Klaus Schulten, and Benoît Roux* J. Chem. Theory Comput. 2017, 13 (12), 5933−5944. DOI: 10.1021/acs.jctc.7b00875

J. Chem. Theory Comput. Downloaded from pubs.acs.org by 185.251.71.111 on 12/08/18. For personal use only.

I

procedure based on the inherent pKa values was to pick a residue according to a fixed probability mass function and then accept or reject a proposal for that residue using eqs 17−19. This procedure was normally repeated until either a new state for any residue was accepted or else a maximum number of (possibly repeated) residues were proposed and rejected. Obviously, the oversight here is that configurations occurring as the result of a rejected change of state should not be discounted and repeating an attempt with a different outcome on a particular residue skews the statistics. The numerical impact of the error was obscured by WHAM, which presupposes that detailed balance is abided and thus imposes a sigmoidal shape. We have rerun our constant-pH MD simulations for dipeptides and obtained slightly different free energy offsets which are needed for simulations of larger systems. These simulations were performed using the new version of the code including corrections for all of the above errors and only moderately different settings (noted below). Most significantly, here we only make one proposal per MC cycle and then decrease the length of the MD portion between cycles in order to increase the number of accepted switches. In the dipeptide simulations, we lower the time between cycles from 10 to 5 ps and in the SNase simulations we lower the time from 10 to 1 ps. For expediency and convenient mapping to computational resources, the SNase simulations were also rerun for 5 ns with eight copies at each pH value. As can be seen in Table 1, with the exception of lysine, all of the new values are within twice the reported standard errors of

n this erratum we wish to report and document a few minor errors in the published algorithm. These errors were repaired soon after their discovery and correct scripts and code were made publicly available (https://www.ks.uiuc.edu/ Research/namd/). The analysis shown here demonstrates that the impact on the published results is minor (at or below statistical precision). The principal error in the original paper concerns eq 19 defining the probability of microscopic protonation states, p(λ), which was subsequently used in a Metropolized independence sampling algorithm (eqs 17 and 18) to propose state-to-state transitions of selected residues based on preassigned inherent pKa values. As written, the text leads one to believe that eq 19 is expressed in terms of the macroscopic inherent pKa value, pK(i),M . However, eq 19 is a strictly correct only if expressed in terms of microscopic pKa values, pK(i),μ a . The difference is that the latter accounts for the multiplicity, ω(λ), of equivalent protonation states for the occupancy vector λ. This multiplicity is just the number of permutations of λ that describe equivalent states and can be written in terms of the number of protons, p, and equivalent protonation sites, q, such that ω = q!/(q − p)!p!. For example, for a protonated carboxylate, ω = 2 (p = 1, q = 2) and for a neutral primary amine ω = 3 (p = 2, q = 3). It follows that (i),M pK(i),μ (λ, λ′) − log ω(λ′)/ω(λ). a (λ, λ′) = pKa Equation 19 can now be more clearly written in terms of microscopic pKa as (i), μ

p(λ) ≡

10 pKa

(λ ′ , λ) − (nλ − nλ′)pH

(i), μ

∑λ ″ 10 pKa

(λ ′ , λ ″) − (nλ ″ − nλ′)pH

(1)

Table 1. Reference Free Energies Computed with the New Code

or alternatively in terms of macroscopic pKa as (i),M

p(λ) ≡

ω(λ)−110 pKa

(i),M

∑λ ″ ω(λ)−110 pKa

ΔF

(λ ′ , λ) − (nλ − nλ′)pH (λ ′ , λ ″) − (nλ ″ − nλ′)pH

(2)

ASP GLU HIS-δ HIS-ϵ CYS LYS

The normalized probability in both of these equations is independent of the reference state λ′ (for clarity, this is no longer an argument to p(λ) as it was in the original paper). In the original paper, ω(λ) was not used explicitly, and the logratio of ω(λ)/ω(λ′) for two states where one proton is transferred was instead written as log[p/(q − p + 1)] (see eq 21 in the original paper). The numerical impact of the error in eq 19 on the simulation results went unnoticed because it was partly canceled out by a programing error. A separate error pertains to our procedure for selecting which residue we should attempt to protonate/deprotonate at each MC cycle. As originally described in the final paragraph under “neMD/MC Sampling” in the original paper, the © 2017 American Chemical Society

old

new

−50.0 −64.4 −0.1 −15.9 −80.8 41.9

−50.3 −64.7 0.1 −15.7 −81.0 42.7

the old values (note that the 95% confidence intervals in the published work were ∼0.4 kcal/mol). The values in most cases deviate only modestly from those previously reported. A similar trend holds for the SNase simulations (Table 2), with the exception of K24, which we had previously noted as a

A

DOI: 10.1021/acs.jctc.8b01075 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Erratum

Journal of Chemical Theory and Computation Table 2. Fitted Apparent pKa Values Reported in the Original Paper (“Old”) and Those Obtained with the New Code Reported Here residue

old

new

GLU

10 43 52 57 67 73 75 101 122 129 135

3.23 4.44 5.01 4.85 4.23 3.48 2.98 4.55 3.90 5.08 3.35

(0.60) (0.07) (0.26) (0.33) (0.80) (0.92) (1.31) (0.45) (0.64) (0.61) (0.48)

3.38 4.68 5.23 4.32 3.78 3.68 2.75 4.36 3.87 4.32 2.91

(0.32) (0.29) (0.33) (0.19) (0.41) (0.31) (0.20) (0.36) (0.51) (0.58) (0.19)

ASP

19 21 40 77 83 95 143 146

2.77 6.78 3.32 0.82 1.97 2.74 4.41 4.01

(0.76) (0.99) (0.52) (0.50) (0.72) (0.39) (0.64) (0.34)

2.75 6.91 3.05 1.81 1.81 2.55 4.56 3.83

(0.24) (0.41) (0.29) (0.59) (0.24) (0.36) (0.19) (0.26)

LYS

24

8.43 (0.45)

>7.2

HIS

8 121

6.66 (0.56) 5.36 (0.50)

6.54 (0.27) 5.13 (0.16)

probably anomalous result. Lastly, the pKa values reported here frequently have smaller standard errors than those previously reported.

B

DOI: 10.1021/acs.jctc.8b01075 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX