Improved Accuracy for Constant pH-REMD Simulations through

Subscriber access provided by Caltech Library

Article

Improved accuracy for constant pH-REMD simulations through modification of carboxylate effective radii Andrew Vincent Yeager, Jason M Swails, and Bill Ray Miller III J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00638 • Publication Date (Web): 14 Sep 2017 Downloaded from http://pubs.acs.org on September 17, 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.

Journal of Chemical Theory and Computation 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 36

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

Journal of Chemical Theory and Computation

Improved Accuracy for Constant pH-REMD Simulations through Modification of Carboxylate Effective Radii Andrew V. Yeager,1 Jason M. Swails,2 Bill Miller III1* 1

Department of Chemistry, Truman State University, 100 E. Normal Ave. Kirksville, MO

63501 2

Department of Chemistry and Chemical Biology and BioMaPS Institute, Rutgers

University, Piscataway, NJ Corresponding author e-mail: [email protected]

Abstract The accuracy of computational models for simulating biomolecules in specific solution pH conditions is critical for properly representing the effect of pH in biological processes. Constant pH (CpH) simulations involving implicit solvent using the AMBER software often incorrectly estimate pKa values of aspartate and glutamate residues due to large effective radii stemming from the presence of dummy protons. These inaccuracies stem from problems in the sampled ensembles of titratable residues that can influence other observable pH-dependent behavior, such as conformational change. We investigate new radii assignments for atoms in titratable residues with carboxylate groups to mitigate the systematic overestimation in the current method. We find that decreased carboxylate radii correspond with increased agreement with experimentally derived pKa values for

ACS Paragon Plus Environment


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

residues in hen egg-white lysozyme (HEWL) and ∆+PHS variants of staphylococcal nuclease (SNase) and improved conformation state sampling compared to experimentally described expectations of native-like structure. Our CpH simulations suggest that decreasing the effective radii of these carboxylate groups is essential to eliminating a significant source of systematic error that hurts the accuracy of both conformational and protonation state sampling with implicit solvent.


Page 2 of 36

Page 3 of 36

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 The pH of a solution can dramatically affect the function of biomolecules. Solution pH may affect protein folding and structure by changing charge distribution within proteins, occasionally destabilizing or denaturing proteins in extreme pH conditions.1,2 Computational models used to simulate biomolecules in specific pH conditions must therefore account for conformational and corresponding behavioral changes that occur. Furthermore, predictions derived from these simulations must accurately reflect experimental findings for simulation to remain a viable mode of study. Traditionally, pH simulations have relied on a fixed protonation state for each titratable residue during simulation, but fixing these protonation states requires prior knowledge of the microscopic pKa of individual residues, which may not be known or measured experimentally. Furthermore, this approach fails to reflect the dynamic nature of protonation states in conditions where multiple protonation states exist in nonnegligible proportions in the desired ensemble.3 When studying the effects of solution pH on a biomolecule, one possible option to overcome the limitations imposed by a modified protonation state model is to allow protonation states to change throughout a molecular dynamics (MD) simulation. Many computational models capable of simulating protonation state changes fall into the category of constant pH molecular dynamics (CpHMD), a set of related models designed to propagate atomic positions through time while simultaneously sampling protonation state changes of titratable residues.4–28 Most popular CpHMD models utilize either a fictitious, continuous titration coordinate for each available state6,8–11,13-17 or utilize Monte Carlo techniques to sample among discrete protonation states.5,7,12,18 Continuous protonation state models track



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

protonation states by adding 'titration particles' to titratable residues sensitive to the solution pH.6,8,14-16 These titration particles, which interpolate between the available protonation states of titratable residues, measure the degree of protonation at those sites on a continuous spectrum ranging from fully deprotonated to fully protonated.6 Discrete protonation state models use metropolized Monte Carlo (MC) exchange attempts at regular intervals during a simulation to sample protonation state changes.5,7,12 The CpHMD model used by AMBER 1419 and in this study relies on a discrete protonation state model that represents protonation state by manipulating partial charges of protonable atoms on titratable residues throughout the MD simulation.7,20 A biomolecule may be simulated using CpHMD to estimate pKa values for individual titratable residues, but such simulations require significant time to converge on a pKa value for a residue.12,20 One option to increase sampling efficiency is to combine CpHMD with replica exchange molecular dynamics (REMD). pH-REMD simulations utilize multiple independent replicas each with a different solution pH that occasionally attempt to exchange pH values.12,21-23 Consequently, coupling REMD with CpHMD enhances protonation state sampling by sampling not only coordinates in a single simulation, but also protonation state space along a pH gradient.18,21,28 Replica exchange operates through periodic MC exchange attempts between independent replicas, allowing adjacent replicas on the pH ladder to exchange when conditions permit.20,28 The probability of a successful replica exchange attempt, based on Eqn. 1, depends on the differences in pH and number of titrating protons between the replicas, where Ni represents the total number of titratable protons in state i.28 ܲ௜→௝ = min{1, ݁‫݈݊[݌ݔ‬10(ܰ௜ − ܰ௝ )(‫ܪ݌‬௜ − ‫ܪ݌‬௝ )]}


(1)

Page 4 of 36

Page 5 of 36

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


If an exchange is successful, simulations continue with the proposed pH until the next replica exchange attempt, while unsuccessful exchange attempts revert to previous conditions, without a replica exchange at that step.18,20,28 The use of either implicit4,7 or explicit9–12 solvent with pH-REMD is supported within AMBER. Explicit solvent simulations result in more reliable pKa predictions in long simulations than implicit solvent simulations.12 Time or computer resource limitations often promote the use of the less computationally expensive implicit solvent models, however, despite the associated reduction in reliability in long simulations. Problems inherent in the model compound definition and underlying force field used in Generalized Born (GB) approximations for implicit solvent simulations can force predicted pKa values away from expected experimental values.12,20 Problems in the underlying force field can lead to unrealistic shapes in the phase space, resulting in inappropriate thermodynamic properties in simulation.7-12,18-20,28-33 Simulations run for no longer than 10 to 20 ns typically yield more accurate pKa predictions than those run up to one or two orders of magnitude longer in implicit solvent.7,12,28,34–36 Unrealistic conformations introduced in longer simulations may change the degree of solvent exposure of titratable residues, degrading pKa estimation over the course of simulations. Implicit solvent simulations also experience faster conformational changes than explicit solvent simulations since the solvent model lacks the viscous effects that slow down conformational sampling in simulations with explicit water, thereby reducing required simulation time.37 Previous work by Swails, et al12 revealed a flaw in the original implementation of pH-REMD simulations stemming from the discrete protonation state model used in the



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

AMBER software suite. The flaw arises from a significant overestimation of the effective radii on carboxylate oxygen atoms caused by the addition of four dummy protons in the syn- and antipositions for both oxygen atoms, shown in Figure 1.7 Since the model treats differences in protonation states as strictly differences in the partial charges of the titratable side-chain, the same effective radii are used for energy calculations involved in the attempted protonation state change. With all four hydrogen atoms included in the calculation of effective radii, the effective radii of the oxygen atoms are systematically overestimated. Swails, et al found that fixing this overestimation is critical to utilizing a GB model for the protonation state changes for explicit solvent CpHMD and pH-REMD simulations, which rely on a hybrid GB-explicit model for MD.12 By reducing the van der Waals radius of oxygen atoms in carboxylate residues from approximately 1.4 Å (original) to 1.3 Å (new), Swails, et al reduced error in pKa estimation of individual carboxylate residues in hen egg-white lysozyme (HEWL) when compared to the original radii.

Figure 1: Positions of all four dummy protons (circled) on carboxylate oxygen atoms of glutamate used in constant pH simulations in AMBER. Similar positions are adopted for aspartate.

The importance of computing the effective radii accurately for GB models has been known for over a decade.38 The GB model used in this study and most frequently in


Page 6 of 36

Page 7 of 36

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


AMBER pH-REMD simulations, developed by Onufriev, Bashford, and Case,39 scales effective Born radii, most significantly in the internal portions of the biomolecule, to control problems that arise from remaining space between atoms within the interior of the molecule.40 Onufriev, Case, and Bashford demonstrated the importance of accurately computing the effective Born radii—with larger values indicating an atom buried more deeply in the protein interior—by calculating the “exact” effective Born radii using the Poisson equation rather than the more commonly employed approximations like the Coulomb field approximation.38 Generally, using the Poisson equation, or the PoissonBoltzmann equation in impure solvent, is impractical due to the high computational cost and general lack of an analytical derivative near the dielectric boundary, which makes GB a far more popular choice for most MD simulations that use implicit solvent.38-39,41 Energy calculations completed using a GB model with ‘perfect’ radii performed better than calculations completed using the effective Born radii estimated by the GB model, suggesting that changes in the computed effective radii can have major effects on the resulting simulation and stressing the importance of accurate radii estimates in GB models.38 Using properly adjusted effective Born radii may therefore help decrease systematic error inherent in the approximations used in GB models, producing results that more closely predict experiment through better conformational sampling and calculation of values derived from sampled ensembles. Nguyen, et al. demonstrated that salt bridge overstabilization involving aspartate residues could be ameliorated by reducing the intrinsic radii from 1.5 Å to 1.4 Å, thereby improving agreement with explicit solvent with their proposed GB-neck2 implicit solvent model.42 The overstabilization of salt bridges is a well-documented problem in most GB



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,43 so it is reasonable to expect reducing the intrinsic solvent radius of the carboxylate oxygen atoms in aspartate and glutamate residues to improve the current GB model as well. However, unlike the typical aspartate and glutamate residues that have no hydrogen atoms bonded to the carboxylate oxygens, the titratable residues defined within AMBER for CpH simulations retain four dummy hydrogen atoms in the syn- and antipositions of both carboxylate oxygens (Fig. 1). Because these atoms are not excluded from the calculation of the effective GB radii—they are excluded only from the chargecharge interactions when in the ‘deprotonated’ state by having their charges set to 0—the carboxylate oxygens feel their presence by introducing a systematic shift in the computed effective radii of the carboxylate oxygens. An increase in the effective radius of the carboxylate oxygen atoms increases the influence of the charges on those atoms by decreasing its dielectric shielding,38 making simulations susceptible to the influence of overall positive or negative charge of a protein. Using the logic of carboxylate radii modification introduced by Swails, et al, this study focuses on the use and performance of GB implicit solvent models in pH-REMD simulations when the same reduction in carboxylate oxygen radii is implemented to compensate for the systematic shift in effective radii. When studying the effects of changing inherent assumptions in pH simulations, hen egg-white lysozyme (HEWL) (Fig. 2) is ideal for testing CpHMD models because it has been extensively studied both theoretically and experimentally and has many titratable residues, including several with a significant shift from “typical” pKa values.2,7,9,12,20,28,44-46 Experimental pKa values for individual residues are available,54 allowing for easy comparison to pKa values predicted by pH-REMD simulations.


Page 8 of 36

Page 9 of 36

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 2: A) Structure of HEWL highlighting the titratable residues shown as sticks: glutamate (red), aspartate (blue), and histidine (purple).B) Structure of ∆+PHS SNase with titratable residues (same colors) and position of mutations (green).

The protein Staphylococcal nuclease (SNase) (Fig. 2) provides a valuable second test as experimental and theoretical pKa estimates of residues within the highly stable ∆+PHS variant are available and experimental pKa estimates indicate that some residues are greatly shifted from standard pKa values.13-15,47-50 Additionally, previous mutation studies detailing the effects of mutant charged residues at a buried position on the stability of the protein across a range of pH values provide an opportunity to explore the ability of CpH systems to predict the conformational consequences of such mutations.47 In this study, we characterize a flaw in the original implementation of the CpHMD model in AMBER affecting carboxylate titratable groups. We investigate the same solution adopted by Swails, et al.12 in the context of implicit solvent simulations to determine whether a similar improvement is seen in the pKa predictions for HEWL . We compare simulations over long time scales for two systems—the first using the original radii



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

assignments for carboxylate oxygens and the second using the modified radii adopted in reference 12. Comparison of implicit and explicit solvent performance for HEWL is completed using the data provided in Swails et al for explicit solvent simulations.12 We further apply the method of Swails, et al12 to three single position mutants and a nonmutant of the ∆+PHS variant of SNase in implicit and solvent to investigate the effect of subtle changes in protein structure on pKa estimates and conformational sampling, simulating both original and modified radii assignments for each system. Short explicit solvent simulations are also completed as described in reference 12 on ∆+PHS SNase variants for structural comparison. We find that utilizing the modified solvation radii decreased error in pKa estimation in most tests and altered structural sampling of the implicit solvent simulations to more closely resemble expectations from experimental data, suggesting that the apparent success of the original radii in calculating the pKas of some carboxylate residues in HEWL and SNase was a fortuitous consequence of error cancellation at the studied time scales. Methods All simulations were run using the HEWL X-ray crystal structure (PDB code 1AKI)51 and mutants of a ∆+PHS SNase variant X-ray crystal structure (PDB code 5KYL).52 The SNase crystal structure, which already contained a deletion of residues 4449 and the P117G, H124L, S128A mutations typical of the stable ∆+PHS variant, was further modified through in silico mutations (G50F, V51N, and V66D, E, and K) to match structures used by Karp, Stahley, and Garcia-Moreno in their study of chargedresidue mutants (V66D, V66E, and V66K).47 The program pdb4amber was used to prepare the PDB files for use with AMBER by changing the names of histidine, aspartate,


Page 10 of 36

Page 11 of 36

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


and glutamate residues to HIP, AS4, and GL4 to match AMBER recognized names for constant pH simulations; add disulfide bonds; and remove crystallographic waters. For HEWL, both original and modified radii implicit solvation pH-REMD simulations were run for a total of 600 ns each, or 50 ns per replicate. For SNase, twelve implicit pHREMD simulations and four explicit pH-REMD simulations were run. Implicit solvent simulations for both original and modified radii were prepared for each of the three mutants (V66D, V66E, and V66K) and the unmutated system – accounting for eight of the implicit systems – and each was run for 600 ns total, or 50 ns per replicate. For the V66K mutant, the mutated lysine was treated as titratable to reflect an experimentallydetermined reduced pKa, though all other lysine residues across simulations were left with fixed protonation states. The remaining four implicit solvent simulations manipulated the identity of titratable residues, with AS4 and HIP residues titratable in two and GL4 and HIP residues titratable in the others. Original and modified radii systems were simulated for each set of titratable residues. Residue 66 was left unmutated for these simulations, which were run for a total of 600 ns each, or 50 ns per replicate. The explicit solvent simulations were completed using the three mutants and the unmutated ∆+PHS variants following the methods described by Swails et al12 for the modification of carboxylate radii in explicit solvent pH-REMD simulations, and run for a total of 120 ns each, or 10 ns per replicate. It is important to note that the reference energy is computed for a particular choice of reference compound parameters. Therefore, changing the starting radii of the carboxylate oxygen atoms requires the reference energies to be recalculated for the AS4 and GL4 residues using the modified solvation radii. These reference energies were



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

defined so that the computed pKa of the reference compound corresponded to its experimental pKa and are summarized in Table 1. Table 1: Reference energies for each protonation state for the titratable carboxylate residues aspartate (AS4) and glutamate (GL4)

Residue AS4 GL4

Original oxygen radii (1.4 Å) reference energies (kcal/mol) Protonated Deprotonated 26.8894581 0 8.4057785 0

Modified oxygen radii (1.3 Å) reference energies (kcal/mol) Protonated Deprotonated 33.299396994 0 15.18076956 0

Each system was run through tleap using the AMBER ff10 underlying force field, which is equivalent to the AMBER ff99SB force field for proteins. The “mbondi2” implicit solvent radii were assigned to each atom. After implicit solvent files were created, the side chain oxygen radii in carboxylate residues (glutamate and aspartate) were changed from the default 1.4 Å (original) to 1.3 Å (new) in the modified radii systems using the parmed utility in AMBER. Using the cpinutil.py utility in AMBER, titratable residues in the original and new radii systems were set to default protonation states, leaving aspartate and glutamate residues deprotonated and the histidine residue protonated. All implicit solvent systems were minimized using 100 steps of the steepest descent algorithm followed by 400 steps of conjugant gradient with a 10 kcal/mol·Å2 restraint force constant applied to the backbone atoms with protonation states held constant. Each system was then heated at constant volume, varying temperature linearly from 10 to 300 K over 1 ns using Langevin dynamics with an atomic collision frequency of 1.0 ps-1. Again, protonation state changes were not attempted. Equilibration extended over 2 ns at constant temperature using Langevin dynamics with an atomic collision frequency of 10.0 ps-1, also at a fixed set of protonation states.


Page 12 of 36

Page 13 of 36

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 explicit solvent systems were minimized using 1000 steps of the steepest descent algorithm followed by 4000 steps of conjugant gradient with a 10 kcal/mol·Å2 restraint force constant applied to water molecules and backbone atoms with protonation states held constant. The systems were then heated from 10 to 300 K at constant volume over 400 ps using Langevin dynamics with an atomic collision frequency of 5.0 ps-1. No protonation state changes were attempted. Equilibration extended over 2 ns at constant temperature using Langevin dynamics with an atomic collision frequency of 1.0 ps-1 with no protonation state changes. All molecular dynamics were run with a 2-fs time step, and a salt concentration of 0.1 M was specified for each GB solvent model, while ions were added in the solvated systems to neutralize charge. Bonds with hydrogen atoms were constrained using SHAKE.41,53 All simulations were run on GPUs using pmemd.cuda32 in AMBER 14.19 The use of Graphics Processing Units (GPUs) to execute pH-REMD simulations allowed for faster MD simulation rates compared to traditional simulations run on Central Processing Units (CPU)s.54-56 In implicit solvent simulations, replica exchange attempts were made every 200 fs, and pH-REMD simulations were run a total of 50 ns per replica. Langevin dynamics were used with an atomic collision frequency of 10.0 ps-1. Protonation state changes were attempted every 10 fs. In the explicit solvent simulations, replica exchange attempts were made every 200 fs with each replica run for a total of 10 ns. Protonation state changes were attempted every 200 fs, with 200 fs of relaxation dynamics. Langevin dynamics were used with an atomic collision frequency of 5.0 ps-1. In all systems, 12 replicas were used in each set of pH-REMD simulations with pH index values ranging from 2.0 to 7.5 and pH interval of 0.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

Simulation analysis required the use of the program cphstats to parse the output files of replicas and obtain protonation state records and calculate predicted pKa values for each residue. The pKa and Hill coefficient (n) were calculated by fitting the deprotonation fraction (fd) and pH for each replica to the Hill equation (Eqn. 2). ଵ

݂ௗ = ଵାଵ଴೙(೛಼ೌష೛ಹ)

(2)

The CpHMD trajectory files were parsed using AMBER’s cpptraj57 analysis package to reorganize trajectory files into ensembles whose snapshots corresponded to a single solution pH to allow for visualization. Root mean square deviations (RMSDs) were also calculated using cpptraj. Comparison between B-factor values and calculated RMSf for α-carbons was also performed for the unmodified HEWL structure. Results Short Peptide Modification To demonstrate the systematic shift in computed effective radii introduced by the presence of four dummy protons on carboxylate residues in pH-REMD simulation, the effective radii calculated for the “normal” aspartate and aspartic acid residues (named ASP and ASH in AMBER, respectively) was compared to that for the titratable aspartate residue (named AS4 in AMBER due to the presence of 4 hydrogen atoms) in Figure 3.


Page 14 of 36

Page 15 of 36

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 3: Effective GB radius is increased by the inclusion of 4 dummy protons on the titratable aspartate residue. Decreasing the intrinsic radius of carboxylate oxygen atoms in the titratable aspartate residue results in an effective GB radius intermediate to the deprotonated aspartate and aspartic acid. Distributions computed over both hydrogen atoms in the model compound sequences ACE-XYZ-NME, where XYZ was either ASH, ASP, or AS4.

Reducing the intrinsic radius of the carboxylate oxygen atoms from 1.4 Å to 1.3 Å results in a distribution of effective radii that is roughly the average of the effective radii of the non-titratable variants, making it an ideal compromise for the AS4 residue. The effect for glutamate carboxylate oxygen is nearly identical to aspartate, and so is not shown. The same arguments apply equally well for reducing the effective radii of glutamate. A more accurate model would account for differences in the effective radii between the different protonation states. However, as this would both hurt computational performance (by requiring additional effective radii calculations for each proposed state) and would require substantial changes to the underlying code, we chose to reduce the intrinsic radii of the carboxylate oxygens by 0.1 Å for this study. Swails, et al. utilized the same reduction in carboxylate oxygen intrinsic radii in the hybrid implicit-explicit solvent model, to great effect, achieving substantial



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 36

improvements in simulation stability and pKa prediction over implicit solvent.12 Thus, any discrepancies in the results between the hybrid approach and this study, then, are almost certainly due to differences in conformational sampling between implicit and explicit solvent models. HEWL A key capability required of any accurate constant pH molecular dynamics model is its ability to predict pKa values for its various titratable residues. Indeed, experimental pKa values are direct reflections of protonation state sampling as a function of solution pH. However, many factors in simulation can affect predicted pKa values including inaccuracies in the reference energies for titratable residues, inaccuracies in the underlying force field leading the simulation to sample unphysical configurations, and inaccuracies in the underlying solvent model which can introduce systematic bias in the protonation state sampling. The predicted pKa values in Table 2 show substantial differences for the two sets of simulations for HEWL (titration plots for implicit simulations available in supporting information (SI): Fig S1-S2).. Table 2: Titratable residue pKa predictions over the first, last, and all 50 ns of simulation (± S.D.) for both the original and new carboxylate oxygen solvent radius, predicted pKas from Swails et al (Expl), and experimental (Expr) pKa values.

Residue Glu 7 His 15 Asp 18 Glu 35 Asp 48 Asp 52 Asp 66 Asp 87 Asp 101 Asp 119 RMSE

Original Carboxylate Radii st 1 ns Last ns 50 ns 3.37 3.90 3.85 ± 0.19 7.26 5.69 5.90 ± 0.64 1.76 2.05 2.09 ± 0.24 5.18 4.78 5.00 ± 0.23 1.55 2.76 2.30 ± 0.48 1.68 2.82 2.47 ± 0.38 1.57 2.41 2.24 ± 0.64 1.75 2.49 2.76 ± 0.30 4.05 2.51 2.54 ± 0.99 2.23 2.04 2.81 ± 0.29 1.07 1.19 1.06

New Carboxylate Radii st 1 ns Last ns 50 ns 3.95 4.30 4.27 ± 0.10 6.49 5.22 5.46 ± 0.24 2.83 3.14 2.92 ± 0.11 6.14 6.22 6.08 ± 0.48 3.73 3.88 3.76 ± 0.06 4.08 4.56 4.58 ± 0.18 4.32 3.63 3.73 ± 0.13 3.63 4.31 4.27 ± 0.15 4.12 3.83 3.92 ± 0.09 3.32 3.58 3.43 ± 0.13 1.43 1.45 1.43


Expl12 120 ns 3.31 6.32 2.89 6.32 1.92 2.63 1.81 2.06 4.25 1.53 0.82

Expr45 2.6 5.5 2.8 6.1 1.4 3.6 1.2 2.2 4.5 3.5 —

Page 17 of 36

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 predicted pKa values appear to degrade with the use of the new radii by comparing the root mean square error (RMSE) only. However, this is misleading as the main contributors to the RMSE for the modified radii come from glutamate 7 and aspartate residues 48, 66, and 87. These four residues are particularly problematic because they are located very near or within highly flexible loop regions of the protein. While explicit solvent treatment reproduces the downward pKa shift of these residues using the same GB model as this study,12 the implicit solvent model drives these residues to increase their solvent accessibility. As a result, their predicted pKas are very similar to those for the model compound. In this case, the GB model’s impact on the protonation state equilibria is expressed through its effect on conformational preferences rather than its direct effect on the electrostatic contributions to the protonation state sampling directly. Residues that are located within more highly structured regions of HEWL, or whose solvent accessibility gives them a small experimental pKa shift—namely His 15, Asp 18, Glu 35, and Asp 52, 101, and 119—show much better agreement with experiment in the modified radii system. In fact, the active site residues Glu 35 and Asp 52 show remarkably better agreement with the modified radii than with the original set. Furthermore, these predictions stay reasonably stable over the course of the full 50 ns simulation, whereas the original radii result in a systematic decrease in the pKa of the general acid (Glu 35) and increase in the pKa of the general base (Asp 52). This effect is caused by the collapse of the active site as described by Swails, et al.12 Comparing the predicted pKas of the two implicit solvent simulations shows that the original radii systematically underestimate the predicted pKas of aspartate and



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

glutamate residues, consistently yielding predicted pKas 1.5 to 2 pK units lower than those for the new radii. This trend was also present when using an explicit solvent model for the MD. However, GB stabilizes unphysical conformational distortions that increase solvent accessibility of each of the carboxylate groups in HEWL. This causes a reduction of the pKa shift from the reference compound, resulting in a systematic increase in predicted pKa for all aspartate residues (except 101 and 119) as well as glutamate 7 and a systematic decrease in the predicted pKa for glutamate 35 toward its reference pKa. When simulated for a long enough period of time, the predicted pKa values using the original radii continue to degrade, as seen both in this study and Swails, et al.12 For short simulations, the conformational distortions do not progress enough to be alarming, and exert just enough effect to result in a fortuitous cancellation of errors in predicting pKas for carboxylate residues in HEWL. By correcting the systematic bias in the electrostatic contributions to protonation state sampling within the GB model, the pKa predictions of well-structured residues improves dramatically at the cost of eliminating the fortuitous error cancellation for some of the less-structured carboxylate residues. The most marked improvement in the predicted pKa comes from the catalytic residue glutamate 35. The prediction with the new radii is remarkably accurate and stable throughout the simulation compared with the substantial underestimate using the original radii that gets progressively worse during the course of 50 ns. The reason for this, as suggested in Reference 12, comes from the collapse of the active site—especially at higher pH values—after long periods of time. In particular, the active site carboxylate residues which serve opposite catalytic function—Glu 35 is the general acid while Asp 52 is the general base67-67—begin to interact very tightly. This effect is demonstrated by


Page 18 of 36

Page 19 of 36

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


plotting the distributions of the Glu 35 – Asp 52 carboxylate group distances as measured from the center of mass of each carboxylate group with each other, shown in Figure 4.

Figure 4: The distribution of Glu 35 – Asp 52 distances is similar throughout the modified radii simulation, but a bimodal distribution is present near the end of the original radii simulation (red dashes), indicating active site collapse.

The initial distributions at the beginning of each simulation appear largely similar between both simulations, with only a very slight preference for smaller distances with the original radii. However, the final 4 ns of each simulation show a very different picture. There are clearly two peaks to the distribution in the simulation with the original radii—one in which the carboxylate groups are separated by less than 5 Å that is absent in the simulation with the new radii. This tight interaction would almost certainly not be observed in the native structure of HEWL, as it would significantly hinder catalytic efficiency of the currently accepted mechanism.44,58 This offers strong evidence that not only do our proposed changes to the solvation radii of the carboxylate oxygens in AMBER improve the quality of protonation state sampling by countering the systematic



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

bias introduced by the dummy hydrogen atoms, it also demonstrably improves conformational sampling of HEWL using an implicit solvent model. Improved conformational sampling does in part seem correlated with increased deviations from starting structures. RMSD values for backbone atoms across replicates were consistently 1 to 2 Å higher in modified radii systems than original, indicating conformational sampling in new radii systems diverged in part by exploring conformations more dissimilar to pre-MD structures than those sampled in original radii systems (SI: Fig S3). A large portion of the disparity seems to stem from increased mobility of less-structured regions in new radii systems, especially in the turn containing Asp 48. The region is highly mobile in both modified and original radii systems, but Bfactor analysis suggests that the region is far more mobile in new radii systems than original (SI: Fig S4). While both system types occasionally deviate from the pattern suggested by the B-factor information from X-ray crystallography, the peak near Asp 48 represents the largest deviation from crystal structure data unshared by original radii systems, and likely accounts for much of the disparity seen in RMSD values between replicates of the same pH. SNase Single residue mutations can greatly change local environments by introducing or removing charge from typically nonpolar or solvent-exposed regions, and accurate constant pH models should be capable of handling the large pKa shifts expected in such cases without forfeiting accuracy of non-mutant residues. Predicted pKa values for aspartate, glutamate, histidine, and mutated residues (V66D, V66E, and V66K) from the


Page 20 of 36

Page 21 of 36

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


modified radii systems performed better than original radii systems when compared to experimental pKas47-50 as shown in Table 3 (SI: Fig S5-S12). RMSE values indicate a general improvement in new radii systems in implicit solvent, as over the entire 50 ns the modified radii systems provided better estimates of pKa overall. One factor which decreased RMSE values for modified systems relative to original was a lack of highly improbable pKa predictions. In multiple cases, original radii systems predicted either negative pKa values or failed to provide an estimate for pKa, despite experimentally determined pKa values falling within the pH range used for replica exchange. New radii systems provided estimates for all residues across simulations and ubiquitously predicted positive values for pKa. Although original radii estimates may have sometimes fallen closer to experimental pKa estimates, disparities between new and original radii systems in such cases were generally small, between 1 and 1.5 pK units, while disparities between modified and original radii predictions were sometimes far larger in instances where modified radii systems provided better estimates, such as estimates for Asp 21 and 40. Across systems, at sites that were not mutated, the new radii systems provided more consistent estimates of pKa original radii systems, potentially indicating a greater degree of consistency in pKa prediction in modified radii systems.


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

Orig.

5.94 ± 0.12

3.52 ± 1.33

1.78 ± 1.03

2.25 ± 110

-5.82 ± 4.14

3.93 ± 0.25

4.08 ± 0.68

4.31 ± 0.15

3.89 ± 0.69

4.26 ± 0.17

3.73 ± 0.16

3.36 ± 0.70

2.37 ± 2.18

-0.02 ± 6.55

3.77 ± 0.18

3.88 ± 0.35

6.30 ± 0.43

3.23 ± 0.47

3.85 ± 0.41

3.26 ± 0.54

2.73

5.10 ± 0.11

3.63 ± 0.38

8.06 ± 0.15

-24.94 ± 186

2.75 ± 0.08

4.60 ± 0.10

5.86 ±0.12

4.62 ± 0.14

4.53 ± 0.28

4.35 ± 0.31

3.88 ± 0.11

0.64 ± 0.66

0.89 ± 0.57

--

3.65 ± 0.16

5.00 ± 0.17

8.18 ± 0.14

3.13 ± 0.26

3.10 ± 0.26

2.96 ± 0.23

7.53

50 ns

Expl.

10 ns

V66D

ACS Paragon Plus Environment 2.41

2.74 ± 0.13

4.46 ± 0.17 1.58

2.27 ± 1.10

3.63 ± 0.23

8.55 ± 0.13

5.19 ± 0.06

4.16 ± 0.23

0.22 ± 37.4

-1.04 ± 0.92

1.73 ± 0.66

3.87 ± 0.24

4.69 ± 0.09

4.32 ± 0.17

4.02 ± 0.14

6.19 ± 0.45

4.04 ± 0.09

3.66 ± 0.12

3.31 ± 0.16

3.56 ± 0.18

4.98 ± 0.44

4.80 ± 0.25

4.20 ± 0.13

4.11 ± 0.48

5.03 ± 0.09

4.25 ± 0.12 3.79 ± 0.62

5.84 ± 0.08

4.58 ± 0.10

4.01 ± 0.13 4.60 ± 0.17

3.28 ± 0.22

-0.04 ± 2.04

4.03 ± 0.08

2.78 ± 0.46

5.29 ± 0.18

Expl.

4.23 ±0.27

5.17 ± 0.41

4.49 ± 0.24

4.30 ± 0.12

5.67 ± 0.11

New

10 ns

1.99

3.52 ± 0.45

4.01 ± 0.39

3.25 ± 0.40

6.44 ± 0.34

3.89 ± 0.23

3.73 ± 0.19

1.45 ± 2.93

2.78 ± 0.49

3.21 ± 0.46

3.77 ± 0.21

4.16 ± 0.16

4.80 ± 0.77

4.25 ± 0.13

3.88 ± 0.27

4.74 ± 0.54

0.85 ± 1.31

--

--

3.82 ± 0.38

6.04 ± 0.11

Orig.

V66E 50 ns

1.29

4.43 ± 0.12

4.46 ± 0.23

4.03 ± 0.24

6.10 ± 0.49

4.25 ± 0.11

3.67 ± 0.2

3.47 ± 0.36

3.78 ± 0.19

2.63 ± 0.32

4.52 ± 0.26

4.28 ± 0.13

5.42 ± 1.68

4.41 ± 0.09

4.71 ± 0.23

4.18 ± 0.17

4.21 ± 0.25

4.99 ± 0.44

3.96 ± 0.27

4.32 ± 0.21

5.82 ± 0.07

New

1.92

2.53 ± 0.25

1.97 ± 1.01

3.45 ± 0.22

8.76 ± 2.19

4.92 ± 0.29

3.99 ± 0.13

0.80 ± 4.92

1.26 ± 7.35

0.94 ± 0.87

4.20 ± 0.16

4.26 ± 0.09

9.12 ± 0.17

4.96 ± 0.15

5.67 ± 0.11

4.66 ± 0.13

2.92 ± 0.16

1.83 ± 1.23

3.99 ± 0.39

3.33 ± 0.21

5.33 ± 0.04

Expl.

10 ns

1.04

3.63 ± 0.54

3.69 ± 0.42

2.71 ± 0.82

6.23 ± 0.82

3.66 ± 0.92

3.80 ± 0.19

1.81 ± 28.0

2.33 ± 1.16

2.68 ± 0.35

3.62 ± 0.26

3.85 ± 0.18

6.84 ± 1.21

4.24 ± 0.94

3.88 ± 0.26

3.49 ± 0.70

1.19 ± 4.94

--

2.61 ± 0.94

3.94 ± 0.32

5.93 ± 0.12

Orig.

V66K 50 ns

0.96

4.41 ± 0.17

4.41 ± 0.12

4.02 ± 0.12

6.05 ± 0.46

4.03 ± 0.12

3.56 ± 0.14

3.56 ± 0.35

3.63 ± 0.21

4.18 ± 0.35

4.58 ± 0.18

3.96 ± 0.19

10.04 ± 1.0

4.18 ± 0.17

4.53 ± 0.20

4.12 ± 0.15

4.01 ± 0.17

5.59 ± 0.45

3.97 ± 0.35

4.20 ± 0.26

5.64 ± 0.12

New

3.76

3.75

3.89

5.4

3.81

2.16

Improved Accuracy for Constant pH-REMD Simulations through

Recommend Documents