
==== Front
bioRxiv
BIORXIV
bioRxiv
2692-8205
Cold Spring Harbor Laboratory

39282392
10.1101/2024.09.03.611076
preprint
1
Article
Force Field Limitations of All-Atom Continuous Constant pH Molecular Dynamics
http://orcid.org/0000-0001-6111-6857
Peeples Craig A.
http://orcid.org/0000-0001-8395-9353
Liu Ruibin
http://orcid.org/0000-0002-3234-0769
Shen Jana
Department of Pharmaceutical Sciences, University of Maryland School of Pharmacy, Baltimore, MD 21201
jana.shen@rx.umaryland.edu
07 9 2024
2024.09.03.611076https://creativecommons.org/licenses/by-nd/4.0/ This work is licensed under a Creative Commons Attribution-NoDerivatives 4.0 International License, which allows reusers to copy and distribute the material in any medium or format in unadapted form only, and only so long as attribution is given to the creator. The license allows for commercial use.
nihpp-2024.09.03.611076.pdf
All-atom constant pH molecular dynamics simulations offer a powerful tool for understanding pH-mediated and proton-coupled biological processes. As the protonation equilibria of protein sidechains are shifted by electrostatic interactions and desolvation energies, pKa values calculated from the constant pH simulations may be sensitive to the underlying protein force field and water model. Here we investigated the force field dependence of the all-atom particle mesh Ewald (PME) continuous constant pH (PME-CpHMD) simulations of a mini-protein BBL. The replica-exchange titration simulations based on the Amber ff19SB and ff14SB force fields with the respective water models showed significantly overestimated pKa downshifts for a buried histidine (His166) and for two glutamic acids (Glu141 and Glu161) that are involved in salt-bridge interactions. These errors (due to undersolvation of neutral histidines and overstabilization of salt bridges) are consistent with the previously reported pKa’s based on the CHARMM c22/CMAP force field, albeit in larger magnitudes. The pKa calculations also demonstrated that ff19SB with OPC water is significantly more accurate than ff14SB with TIP3P water, and the salt-bridge related pKa downshifts can be partially alleviated by the atom-pair specific Lennard-Jones corrections (NBFIX). Together, these data suggest that the accuracies of the protonation equilibria of proteins from constant pH simulations can significantly benefit from improvements of force fields.
==== Body
pmcIntroduction

Solution pH mediates many important biological processes through coupling proton titration with conformational changes of proteins. Over the last two decades, constant pH methods have been developed to rigorously account for solution pH in condensed-phase molecular dynamics (MD) simulations. Instead of fixing the protonation states of protein sidechains to the initial conditions, e.g., histidines with one proton on either the δ or ϵ imidazole nitrogen, constant pH simulations allow protonation states to respond to changes in the electrostatic environment accompanying conformational dynamics. Currently, the main approaches to enable proton-coupled dynamics is through Monte-Carlo (MC) sampling of protonated and deprotonated states (discrete or hybrid MD/MC constant pH methods)1–5 and an extended Hamiltonian description whereby an auxiliary set of fictitious particles are propagated to represent proton titration (continuous constant pH or λ dynamics based constant pH methods).6–11 A more detailed discussion of the development of constant pH methods is given in a recent review12 for more references.

In the early developments of constant pH methods, various generalized Born (GB)2,6,7,13 or Poisson-Boltzmann (PB) implicit solvent models are used for both conformational and protonation state sampling. The so-called hybrid-solvent scheme combines the conformational sampling in explicit solvent with implicit-solvent calculation of titration energies4,14 or forces,8 resulting in significantly improved accuracies of the calculated pKa values. The more recent developments11,15–18 utilize explicit solvent for both conformational and protonation state sampling, which allows constant pH simulations to study systems that implicit-solvent representations are either insufficiently accurate or simply unfeasible to describe, for example, highly charged systems or those in a heterogeneous dielectric environment. The first all-atom λ dynamics based CPHMDMSλD implementation in CHARMM was applied to calculate the pKa’s of RNAs19 and peptides inserted in the lipid bilayer.20

By removing the deficiencies due to the implicit-solvent models, all-atom constant pH simulations can in principle give more accurate pKa values compared to the hybrid-solvent constant pH simulations. The benchmark simulations using the first all-atom particle mesh Ewald (PME) continuous constant pH MD (CpHMD) implementation in CHARMM16 demonstrated improved correlation between the calculated and experimental pKa shifts relative to the solution (also known as the model) values as compared to its “predecessor”, the hybrid-solvent CpHMD in CHARMM,8 even though the root-mean-square errors (rmse) with respect to the experimental data are similar.

The key physical determinants of protein pKa shifts relative to the model values are electrostatic interactions and desolvation free energies, which have been shown to vary between different protein force fields for proteins and water.21–26 As such, pKa values derived from constant pH simulations can be used to probe the force field deficiencies.27 So far, the benchmark simulations of all-atom constant pH methods15–18,28 have been conducted using the CHARMM c22/CMAP,29,30 c36,31 or c36m32 force fields, while other force fields such as the widely used Amber ff14SB33 or ff19SB34 force fields have not been tested. Recently, Sequeira et al. compared the pKa calculations using the hybrid-solvent MD/MC based constant pH simulations with the GROMOS 54A7 and CHARMM c36m force fields.35 Their study showed that, with the c36m force field the calculated pKa’s are very stable over the simulation time (up to 50 ns), whereas with the 54A7 force field the pKa’s drift significantly over the simulation time.35

In this work, we investigate the force field dependence of the all-atom PME CpHMD simulations by comparing the titration simulations of a mini-protein BBL (Figure 1) based on the Amber ff19SB,34 ff14SB,33 and CHARMM c22/CMAP29,30 protein force fields and their respective preferred water models, TIP3P,36 OPC,37 and CHARMM-style TIP3P.29 BBL is a well suited benchmark protein for the constant pH based pKa calculations, as it contains buried histidines and carboxylic acids involved in salt-bridge interactions, for which large pKa shifts relative to solution (also called model) values are expected and challenging to predict accurately.38 Also, BBL has been previously used to benchmark the performance of the all-atom PME CpHMD implementations in CHARMM39 and Amber programs.40 Importantly, due to its small size, conformational sampling is less likely an issue. We found that the largest error for all force fields is for a buried histidine, while glutamic acids involved in salt bridge interactions have overestimated pKa shifts. Specific ion binding was observed in the ff19SB simulations; however, the use of atom-pair specific corrections to the Lennard Jones parameters (widely known as NBFIX)41–45 did not reduce the errors of the calculated pKa’s.

Results and Discussion

Comparison of the calculated pKa’s based on the Amber and CHARMM force fields and comparison to experiment.

BBL has 6 carboxylic acids and 2 histidines. The calculated pKa’s based on the ff19SB protein force field34 gave a root-mean-square error (rmse) of 1.26 pH units with respect to experiment, which is much larger than the rmse of 0.62 with the CHARMM c22 protein force field29,30 and the rmse of 0.66 from the GB-Neck253 implicit-solvent simulations with the ff14SB force field33 (Table 1). Convergence and titration curves of ff19SB simulations are given in Supplemental Figure S1–S4. Note, the rmse of the pKa’s calculated from the preliminary ff14SB simulations is 1.34, similar to ff19SB (Supplementary Table S1). Curiously, for both ff19SB and c22 simulations, the largest pKa calculation error is for His166, which has the respective calculated pKa’s of 2.4 and 4.1, corresponding to the pKa downshifts of 4.1 and 2.4 relative to the model value of 6.5 (Table 1). Both ff19SB and c22 simulations correctly predicted the direction of the pKa shift; however, the magnitude of the downshift is too large by 3.0 and 1.2 pH units as compared to experiment (Table 1). For the ff14SB force field, the downshift is overestimated by 2.6 unit (Supplementary Table S1). Without His166, the rmse for the ff19SB and c22 simulations are reduced to 0.73 and 0.48, respectively (Table 1), while for ff14Sb the rmse without His166 is 0.99 (Supplementary Table S1). By contrast, in the GBNeck2 implicit-solvent simulations, the pKa downshift of His166 is underestimated by by 0.6.

Following His166, the largest pKa calculation errors with the ff19SB force field are for Glu141 and Glu161, which have the calculated pKa’s of 3.4 and 2.5, corresponding to the pKa downshifts of 0.8 and 1.7 relative to the model pKa of 4.2, respectively (Table 1). With the c22 force field, the calculated pKa’s for Glu141 and Glu161 are 4.0 and 4.0, corresponding to the pKa downshifts of 0.2 relative to the model value (Table 1). Thus, the calculated pKa shifts are in the same direction with the two force fields. Compared to the experimental pKa’s of 4.5 and 3.7 for Glu141 and Glu161, the ff19SB simulations underestimated both pKa’s by about 1.1 unit, whereas the c22 simulations underestimated the pKa’s by 0.5 and 0.3 units, respectively. Below we discuss the origins of the underestimated pKa’s for Glu141, Glu161, and His166.

pKa downshifts of Glu141 and Glu161 are due to salt bridge formation.

The pKa’s of Glu141 and Glu161 calculated by the ff19SB simulations are both downshifted with respect to the model value, suggesting that they are involved in attractive electrostatic interactions. Trajectory analysis shows that Glu141 occasionally forms a salt bridge with Arg137, with the occupancy (or probability) increasing from zero at pH 1 to a maximum of nearly 20% near pH 4 where deprotonation is completed (Figure 2a–c). Over the same pH range, the solvent exposure of Glu141, defined by the number of water in the first solvent shell of the carboxylate oxygens, also increases and plateaus to about 7 upon complete deprotonation (Figure 2d). Since the first solvent shell of carboxylate oxygens of the model pentapeptide (AAEAA) contains an average of 7 water molecules (SI), these analyses (Figure 2c,d) suggests that Glu141 becomes fully exposed to solvent once the salt bridge is disrupted.

The strong correlation between the sigmoidal shaped pH profiles of the deprotonation fraction, salt-bridge occupancy, and solvent exposure of Glu141 indicates that the salt-bridge interaction, which stabilizes the charged state, is the main determinant for the pKa downshift. This phenomenon has been observed in the previous titration simulations using the CpHMD implementations in the CHARMM and AMBER programs with the c22 force field.16,28 However, since the experimental pKa of Glu is 4.5, which is 0.3 up shifted relative to the model value, the Glu141–Arg137 salt bridge interaction is likely overly stabilized by the ff19SB force field, and the overstabilization is to a lesser extent with the c22/CMAP force field as the pKa underestimation is 0.5 unit smaller.

We next examined Glu161, which shows a pKa downshift that is 1.2 units too large compared to experiment (Table 1). Similar to Glu141, deprotonation of Glu161 is also strongly correlated with the salt-bridge formation (with either Lys165 or Arg160, Figure 3a–d). The salt-bridge occupancy of Glu161 increases to a maximum of 28% near pH 5 (Figure 3d), which is higher than Glu141 (Figure 2c). Over the same pH range, the number of water in the first solvent shell of Glu161 increases to about 6.5 (Figure 3e), which is slightly lower than that for Glu141 (Figure 2d). The increased salt-bridge interactions and slightly decreased solvent exposure favor the charged state, which may explain the lower pKa of Glu161 relative to Glu141.

Downshifted pKa of Asp162 is due to hydrogen bonding.

The calculated pKa’s of Asp162 based on the ff19SB and c22 force fields are similar and downshifted relative to the model value by 0.4 and 0.7, respectively. These pKa downshifts are similar to the experimental downshift of 0.5 units. Trajectory analysis suggests that the pKa downshift of Asp162 can be attributed to the formation of a hydrogen bond (h-bond) with Thr152 (Figure 4a). The deprotonation of Asp162 occurs over the pH range of 1 to 5 (Figure 4b), where the increased deprotonation is accompanied by the h-bond formation between the carboxylate of Asp162 and the hydroxyl group of Thr152 which increases from 5.3% to 83.8% (Figure 4c). Stabilization of the deprotonated carboxylate by accepting a h-bond has been frequently observed in the c22 based CpHMD simulations using the CHARMM and Amber programs.16,28 The similarity between the experimental and calculated pKa’s of Asp162 based on the ff19SB and c22 force fields suggests that the h-bond interaction involving Asp162 may be represented appropriately by the force fields.

pKa downshift of His166 is due to solvent sequestration.

To understand the significant pKa downshift of His166, we examined its solvent exposure and hydrogen bond (h-bond) environment in the titration simulations (Figure 5a,b). Consistent with the analysis of the c22 titration simulations,28 the pKa downshift of His166 can be mainly attributed to solvent displacement. Above pH 3.5, the first solvent shell of His166 includes only 3 water, As His166 becomes protonated with decreasing pH, the number of water in the first solvation shell of His166 increases from about 3 at pH 4 to about 5 at pH 1 (Figure 5c,d). Consistent with the c22 titration simulations,28 the h-bond formation of His166 is minimal. At pH 1–2, His166 (at Nδ, Supplemental Figure S5) donates a h-bond to the backbone carbonyl oxygen of Asp162, with an occupancy of ∼11% (Figure 5e).

NBFIX improves the pKa calculations.

To reduce the binding interactions between ions and proteins and between oppositely charged ions, atom-pair specific corrections to the Lennard-Jones parameters (commonly known as NBFIX) have been introduced to both CHARMM41,42 and Amber43–45 force fields. The ff19SB titration simulations showed specific binding of Cl− with arginine and histidine and Na+ with glutamic acid (see below), and such interactions were not present in the c22 simulations28 with the NBFIX corrections.42

Since the NBFIX corrections for ff19SB have not been developed, we repeated the simulations of BBL using the NBFIX corrected ff14SB force field (ff14SBfix)43–45 and compared to the pKa’s calculated in the previous work28 based on the ff14SB force field, which gave a rmse of 1.78 or 1.38 without His166 (Supplemental Table S1). These errors are larger than those based on the ff14SBfix simulations, which gave a rmse of 1.58 or 1.14 without His166 (Supplemental Table S1). It is worth noting that even with the NBFIX corrections, the pKa errors based on ff14SB are still larger than the ff19SB simulations without the NBFIX corrections (rmse of 1.26 or 0.68 without His166). Below we discuss ion binding to the aforementioned residues, Glu141, Glu161, and Asp162, and His166 in the ff19SB simulations and the effects of using the NBFIX corrections.

NBFIX abolishes ion binding and weakens salt bridges involving Glu141 and Glu161.

In the case of Glu141, the ff19SB simulations showed a considerable amount of Cl− binding to Arg137, which decreases from about 30% at pH below 4 to about 19% at pH 7.0 (Figure 2a and 2e). There is also occasional binding of a Na+ ion with Glu141, with a maximum occupancy of about 9% (Figure 2a and 2e). The introduction of NBFIX significantly weakened the interactions between Cl− and Arg137 and between Na+ and Glu141, with the highest occupancy below 10% for either ion binding in the entire simulation pH range (Figure 2e). At the same time, the salt-bridge interaction between Glu141 and Arg137 is nearly abolished (Figure 2c), resulting in the increased solvent accessibility of Glu141 (Figure 2d). Comparing the pKa of Glu141 based on ff14SBfix and ff14SB, the 0.2-unit upshift can be attributed to the weakened salt bridge (Supplementary Table S1).

In the case of Glu161, the ff19SB simulations showed significant occupancies of chloride binding to Arg160 in the entire simulation pH range, which were abolished in the ff14sbfix simulations (Figure 3f). At the same time, the salt-bridge interactions of Glu141 with Lys165 and Arg160 were also disrupted (Figure 3d), while the solvent accessibility was slightly increased (by half of a water, Figure 3e). We suggest that the abolishment of the salt-bridge interactions is the major reason for the 0.4-unit reduction in the pKa downshift for Glu161 in the ff14sbfix as compared to the ff14SB simulations (Supplementary Table S1).

The exaggerated pKa downshift of Asp162 is reduced by the ff19SB relative to the ff14SB force field.

As to Asp162, the ff19SB and ff14SBfix simulations showed no ion binding or salt bridge formation. This explains why the calculated pKa’s based on the ff14SBfix and ff14SB force fields are nearly identical. However, they are respectively 1.6 and 1.7 units lower than the value from the ff19SB simulations (Table 1 and Supplementary Table S1). The lower pKa may be attributed to the stronger h-bond interaction with the hydroxyl group of Thr152, which stabilizes the deprotonated state of Asp162 (Figure 4c). Consequently, the pH profile of the h-bond occupancy is shifted to a lower pH range, coinciding with the pH range of the aspartic acid deprotonation (Figure 4b).

NBFIX removes chloride binding with His166 and backbone h-bonding may be overstabilized by ff14SB.

In the ff19SB simulations, a chloride ion occasionally binds His166 (mostly at the imidazole Nδ), with an occupancy increasing from about 5% above pH 4 to over 15% below pH 2 (Supplementary Figure S5). The pH dependence of ion binding is inversely correlated with the histidine deprotonation, which is consistent with the electrostatic attraction between the positively charged histidine and Cl− at low pH. The ff14SBfix simulations abolished ion binding (Figure 5f), and increased the pKa of His166 by 0.2 units relative to the ff14SB simulations (Supplementary Table S1). Interestingly, the pKa based on ff19SB is 0.3 units higher than ff14SBfix (Table 1). This may be attributed to the increased solvent accessibility, particularly below pH 4 where the number of water is 5 in the ff19SB as compared to 2 in the ff14SBfix simulations; this increased solvation stabilizes the charged state of histidine resulting in a higher pKa Ḟurthermore, the h-bonding interaction between His166 (exclusive to Nδ) and the backbone carbonyl of Asp162 is negligible in the ff19SB simulations, whereas its occupancy increases with pH to nearly 30% at pH 8 in the ff14SBfix simulations (Figure 5e). It is conceivable that the increased solvent exposure is due to the disruption of the sidechain-to-backbone h-bond.

Concluding Discussion

In this work, we investigated the force field dependence of the all-atom PME CpHMD simulations using the pKa calculations for a mini-protein which contains downshifted pKa’s of several carboxylic acids and one histidine. The CHARMM c22 and Amber ff19SB and ff14SB force fields were considered. Our data showed that the ff19SB and ff14SB force fields overestimate the pKa downshifts of Glu141 and Glu161 involved in salt-bridge interactions and the pKa downshift of the buried His166. These trends are consistent with the c22 force field but the overestimation by the ff19SB and ff14SB force fields is to a greater extent.

The pKa data and pH-dependent conformational analysis suggest that the salt-bridge interactions involving Glu may be overly stabilized, which is a known deficiency of protein force fields.23,26,45,54 Application of the NBFIX corrections that significantly weaken ion binding and salt-bridge interactions in the ff14SB simulations43–45 demonstrated small upshifts of the calculated pKa’s; however, the reduction in the pKa downshift remains insufficient. This is a topic that deserves further investigation in the future. Compared to the all-atom simulations, the pKa’s of Glu141 and Glu145 from the GBNeck2 simulations are much closer to experiment, although the pKa downshifts remain slightly overestimated. This improvement may be in part explained by the larger Lys+–Glu− and Arg+–Glu− salt-bridge distances due to the difference in the geometry between GB and TIP3P.53

The overestimation of the pKa downshift of His166 from the ff19SB and ff14SB simulations can be explained by the underestimated hydration free energy of the neutral histidine25 by the ff14SB force field and TIP3P or OPC355 (similar to OPC37) water model, which suggests that the desolvation free energy of a buried histidine is too large. Interestingly, the pKa of His166 is improved based on ff19SB relative to ff14SB, which may be attributed to the weakened intramolecular h-bonding and increased solvent accessibility (stronger solute-water interactions) in the OPC water.37 Compared to the Amber force fields, the hydration free energy of the neutral histidine based on the c22 force field and TIP3P model is closer to experiment, although that of the Hid tautomer is also underestimated.56 The difference in the hydration free energies between the ff14SB and c22 force fields explains why the pKa of His166 is shifted lower based on the ff19SB or ff14SB relative to the c22 simulations. In stark contrast to the all-atom CpHMD simulations, the pKa of His166 is somewhat overestimated in the GBNeck2 simulations, which may be due to the increased solvent exposure as a result of larger conformational fluctuation in the GB-Neck2 solvent.57

The largest improvement between the ff19SB and ff14SB calculated pKa’s is for Asp162. While ff19SB gives a pKa in agreement with experiment, the ff14SB or ff14SBfix overestimates the pKa downshift by 1.5 or 1.6 units. Similar to His166, the h-bond involving Asp162 is weakened in the ff19SB simulations, which may be attributed to solute-water interactions in the OPC water37

Taken together, this study confirms that the accuracy of pKa calculations using constant pH simulations is dependent on the underlying protein force field and water model. Given sufficient sampling, for example, in the case of the mini-protein BBL, deviations between the calculated and experimental pKa’s are largely reflective of the limitations of the force fields/water models. Salt-bridge overstabilization and underestimation of the histidine neutral state hydration are two common deficiencies of the ff19SB, ff14SB and c22 force fields (combined with their respective water models), although the magnitude of such errors appears to be larger for the ff19SB or ff14SB force fields. The BBL data also confirms that ff19SB with OPC is more accurate than ff14SB with TIP3P in terms of the balance between intramolecular and intermolecular protein-water interactions. Our work represents an initial effort to understand the force field limitations with the goal of improving the accuracy of all-atom PME CpHMD simulations.

Methods and Protocols

Unless otherwise noted, the simulations were performed using the all-atom particle-mesh (PME) CpHMD implementation28 in Amber24.40

Preparation of the protein system.

The coordinates of BBL were retrieved from the Protein Data Bank (PDB) entry 1W4H.46 The first NMR model was used. Following our previous protocol,58 the positions of hydrogens were built using the HBUILD command in the CHARMM program,39 and a custom CHARMM script was used to add dummy hydrogens to the syn positions of carboxylate oxygens on all Asp and Glu sidechains. Next, the CHARMM coordinate file was converted to the Amber format. The protein was solvated in a truncated octaheral water box with a minimum of 15 Å between the protein heavy atoms and the water oxygens at the box edges. To neutralize the simulation box at pH 7.5 and reach an ionic strength of 150 mM, 18 Na+ and 19 Cl− were added.

Minimization, heating, and equilibration of the protein.

The pmemd.cuda engine of the Amber2024 program40 was used for simulations. First, the protein system underwent 10,000 steps of energy minimization using the steepest descent (1000 steps) and conjugate gradient (9000 steps) algorithms, whereby the protein heavy atoms were harmonically restrained with a force constant of 100 kcal·mol−1·Å−2. Next, the PME-CpHMD module28 was turned on and the restrained system was heated from 100 to 300 K over 1 ns at pH 7.5 with a 1 fs timestep in the NVT ensemble. Finally, a two-stage equilibration was performed in the NPT ensemble. In the first stage, the PME-CpHMD simulation was performed at pH 7.5 with a 1-fs time step. A harmonic restraint was placed on heavy atoms, and it was reduced from 100 in the first 1 ns to 10 kcal·mol−1·Å−2 in the second 1 ns. In the second stage, 16 pH replicas were placed with increments of 0.5 units in the pH range 1.0–8.5 and simulated independently for 2 ns, employing a 2 fs timestep and harmonic restraints on the backbone heavy atoms. Here the restraints were gradually reduced in four 500-ps steps, from 10 to 5, 2.5, and 0 kcal·mol−1·Å−2.

Production pH replica-exchange CpHMD simulations.

In the production run, the asynchronous pH replica-exchange (REX) protocol59 was turned on using the same 16 pH replicas as in the equilibration stage. The exchange between adjacent pH was attempted every 1000 MD steps (2 ps) according to the Metropolis criterion.8 The REX simulations were run for 32 ns per replica, with the aggregate simulation time of 512 ns. The pKa’s for all residues were converged (Supplementary Figure S1–S4). The first 10 ns per replica was discarded in the conformational analysis and pKa calculations.

MD protocol.

The SHAKE algorithm was employed to constrain the bonds involving hydrogens whenever a 2-fs step was used. The production runs were performed in the NPT ensemble, where the temperature was maintained at 300 K using the Langevin thermostat with a collision of 1.0 ps−1, and the pressure was maintained at 1 atm using the Berendsen barostat with a relaxation time of 1 ps. The PME method was used to calculate long-range electrostatics with a real-space cutoff of 12 Å and a 1-Å grid spacing in the reciprocal space calculation. Consistent with the our previous titration simulations with the c22 force field,28 Lennard-Jones energies and forces were smoothly switched off over the range of 10 to 12 Å.

Model pKa values and pKa calculations.

The model titration parameters for ff19SB and ff14SB simulations of Asp, Glu, and His were taken from the previous work.28 The validation simulations28 were conducted using titration simulations of penta-peptides ACE-AlaAlaXAlaAla-NH2 (X = Asp, Glu, or His) at independent pH conditions with an interval of 0.5 pH in the pH range 2–5.5 for Asp, 2.5–6 for Glu, and 4–8 for His. The simulations lasted 20 ns at each pH. For ff19SB, the calculated pKa’s based on bootstrap from the deprotonated fractions at different pH are 3.32±0.07, 3.92±0.05, and 6.48±0.06 for Asp, Glu, and His, respectively28. For ff14SB, the calculated pKa’s are 3.51±0.11, 4.07±0.07, and 6.71±0.10 for Asp, Glu, and His, respectively.28 Given the deviations from the target NMR derived pKa’s60 of 3.67, 4.25, and 6.54 for Asp, Glu, and His penta-peptides, we made the following corrections to the protein pKa’s. For ff19SB simulations, the corrections were 0.35 for Asp, 0.33 for Glu, and 0.06 for His. For ff14SB simulations, the corrections were 0.3 for Asp and Glu is 0.0. Note, simulations of Asp and Glu penta-peptides with the NBFIX corrections gave no different pKa’s as those from the ff14SB simulations.

We previously showed that protein pKa’s derived from the PME-CpHMD simulations are affected by the simulation box size, and the finite-size effect decreases with increasing number of water, becoming negligible with large water boxes.28 Since the current simulations employed a very large water box (15 Å cushion between protein and box edges), the calculated finite-size16,28 corrections were small (no larger than 0.1 unit). As such, no finite-size corrections were made to the calculated pKa’s.

Force field parameters and CpHMD specific changes.

The peptides and proteins were represented by the Amber ff19SB34 or ff14SB33 protein force field, with the respective OPC37 or TIP3P50 model for representing water. The ion parameters were taken Ref.49 For the simulations with ff14SB, the NBFIX corrections developed by Yoo and Aksimentiev (CUFIX) were applied to the atom pair specific Lennard-Jones parameters to destabilize the attractive interactions between Na+ and carboxylate (Asp− or Glu−),43 between amine (Lys+) or guanidinium (Arg+) and carboxylate (Asp− or Glu−).44,45 The all-atom PME CpHMD implementations in Amber include two force field modifications common to constant pH methods.2,13,28 First, to allow the single reference concept in constant pH simulations, the partial charges on the backbone are fixed to the values of a single protonation state (Asp− /Glu− and Hid/Hie), and the residual charge (ranging from 0.10 to 0.14 e for Asp, Glu and His) is absorbed onto the C-β atom for titration dynamics.13 Second, the rotation barrier around the carboxylate C-O bond is increased to 6 kcal/mol to prevent the dummy proton from rotating from the syn (initial position in the set up) to the anti position, in which case the proton would lose the ability to titrate.2,13

Supplementary Material

Supplement 1

Acknowledgement

Financial support is provided by the National Institutes of Health (1R35GM148261 and R01CA256557).

Data Availability

Parameter and simulation input files are freely available on github.com/JanaShenLab/CpHMD_ff_comparison.

Figure 1: Structure and titratable sites of BBL.

The NMR model of BBL (PDB ID: 1W4H,46 first entry) The titratable sidechains are labeled and shown in the stick model.

Figure 2: Titration of Glu141 in BBL and the pH-dependent salt bridge formation, solvent exposure, and ion binding.

a) Snapshot of the bound ions from the ff19SB trajectory at pH 6.5. Na+ and Cl− ions are shown as purple and yellow spheres, respectively. b) pH-dependent deprotonated fraction of Glu141. c) pH-dependent occupancy of the salt bridge formation between Glu141 and Arg137. A salt bridge is considered formed if the distance between a carboxylate oxygen of the charged Glu141 and a guanidinium nitrogen is below 4.0 Å. d) Number of water molecules within 3.4 Å of a carboxylate oxygen of Glu141 at different pH. e) Occupancy of Cl− binding of Arg137 (closed circles) and Na+ binding of Glu141 (open circles) at different pH. An ion is considered bound if it is within a cutoff distance, 4.2 Å for Cl− and 4.0 Å for Na+, from a guanidinium nitrogen of Arg137 or a carboxylate oxygen of Glu141. The data based on the ff19SB and ff14SBfix force fields are shown in red and black, respectively, and do not include pKa corrections.

Figure 3: Titration of Glu161 in BBL and the pH-dependent salt bridge formation, solvent exposure and ion binding.

a) Snapshot showing a salt bridge between Glu161 and Arg160 taken from the ff19sb trajectory at pH 7. b) Snapshot of Arg160 bound to Cl− and Glu161 forming a salt bridge with Lys165 from the ff19sb trajectory at pH 7. c) Deprotonated fraction of Glu161 at different pH. d) Salt bridge occupancy of Glu161 with either Arg160 or Lys165 as different pH. e) Number of water within 3.4 Å of any carboxylate oxygen of Glu161 at different pH. f) Occupancy of Cl− binding to Arg160 at different pH. The data based on the ff19SB and ff14SBfix force fields are shown in red and black, respectively, and do not include pKa corrections.

Figure 4: Titration of Asp162 in BBL and the pH-dependent hydrogen bonding.

a) Snapshot of a hydrogen bond between Thr152 and Asp162 from the ff14sbfix trajectory at pH 7. b) pH-dependent deprotonation fraction of Asp162. c) pH-dependent occupancy of the hydrogen bond (h-bond) formation between the carboxylate of Asp162 and the hydroxyl group of Thr152. A h-bond is considered formed if the distance between the donor and acceptor heavy atoms is below 3.5 Å and the donor−H···acceptor angle is greater than 150°. d) Number of water within 3.4 Å of a carboxylate oxygen of Asp162 at different pH. The data based on the ff19SB and ff14SBfix force fields are shown in red and black, respectively, and do not include pKa corrections.

Figure 5: Titration of His166 in BBL and the pH-dependent hydrogen bonding, solvent exposure, and ion binding.

a) Hydrogen bond formation between the imidazole nitrogen of His166 and the backbone carbonyl oxygen of Asp162 in a snapshot taken from the ff19SB trajectory at pH 2.5. b) The deprotonation fraction of His166 at different pH. c) Occupancy of the hydrogen bond (h-bond, shown in a) between His166 and Asp162. A h-bond is defined by a maximum distance of 3.5 Å between N and O and a minimum angle N–H O of 150°. d) The imidazole nitrogen on His166 is bound to Cl− (yellow sphere), shown in a snapshot taken from the ff19sb trajectory at pH 2.5. e) The number of water within 3.4 Å of an imidazole nitrogen on His166 at different pH. f) The occupancy of Cl− binding of His166 at different pH. A binding event is defined as a Cl− ion within 4.2 Å of an imidazole nitrogen on His166. The data based on the ff19SB and ff14SBfix force fields are shown in red and black, respectively, and do not include pKa corrections.

Table 1: Comparison of the calculated pKa’s of BBL using the all-atom PME CpHMD titration based on Amber and CHARMM force fieldsa

Residue	Expt	ff19sb	ff14sbfix	c22	ff14sb/GB	
D129	3.9	3.5	3.6	3.6	2.8	
E141	4.5	3.4	3.1	4.0	4.1	
H142	6.5	6.6	6.5	5.8	6.9	
D145	3.7	3.3	2.2	3.2	2.6	
E161	3.7	2.6	2.6	4.0	3.3	
D162	3.2	3.3	1.7	2.9	3.2	
E164	4.5	3.6	3.4	4.0	4.0	
H166	5.4	2.4	2.1	4.2	6.0	
maxe		−3.0	−3.3	−1.2	−1.1	
rmse		1.26	1.58	0.61	0.66	
rmse (w/o H166)		0.73	1.14	0.48	0.67	
a Maximal error (maxe) and root-mean-square error (rmse) are calculated with respect to the experimental data. The experimental data of BBL is taken from Ref. 47,48 ff19sb refers to the ff19SB protein force field 34 with the OPC water model 37 and Amber default ion force field. 49 ff14sbfix refers to the ff14SB protein force field 33 with TIP3P water model 50 and the NBFIX corrected ion force fields. 43–45 c22 refers to the CHARMM c22 protein force field 29 with the CHARMM style TIP3P water model 29,50 and the NBFIX corrected ion force fields. 42,51,52 Column ff14sb/GB contains the previous data 13 obtained from the pH replica-exchange titration simulations based on the ff14SB force field and GB-Neck2 implicit-solvent model. 53 Residues discussed in the main text are highlighted in bold font.
==== Refs
References

(1) Baptista A. M. ; Martel P. J. ; Petersen S. B. Simulation of Protein Conformational Freedom as a Function of pH: Constant-pH Molecular Dynamics Using Implicit Titration. Proteins 1997, 27 , 523–544.9141133
(2) Mongan J. ; Case D. A. ; McCammon J. A. Constant pH Molecular Dynamics in Generalized Born Implicit Solvent. J. Comput. Chem. 2004, 25 , 2038–2048.15481090
(3) Stern H. A. Molecular Simulation with Variable Protonation States at Constant pH. J. Chem. Phys. 2007, 126 , 164112.17477594
(4) Swails J. M. ; York D. M. ; Roitberg A. E. Constant pH Replica Exchange Molecular Dynamics in Explicit Solvent Using Discrete Protonation States: Implementation, Testing, and Validation. J. Chem. Theory Comput. 2014, 10 , 1341–1352.24803862
(5) Chen Y. ; Roux B. Constant-pH Hybrid Nonequilibrium Molecular Dynamics–Monte Carlo Simulation Method. J. Chem. Theory Comput. 2015, 11 , 3919–3931.26300709
(6) Lee M. S. ; Salsbury F. R. ; Brooks C. L. III Constant-pH Molecular Dynamics Using Continuous Titration Coordinates. Proteins 2004, 56 , 738–752.15281127
(7) Khandogin J. ; Brooks C. L. III Constant pH Molecular Dynamics with Proton Tautomerism. Biophys. J. 2005, 89 , 141–157.15863480
(8) Wallace J. A. ; Shen J. K. Continuous Constant pH Molecular Dynamics in Explicit Solvent with pH-Based Replica Exchange. J. Chem. Theory Comput. 2011, 7 , 2617–2629.26606635
(9) Donnini S. ; Tegeler F. ; Groenhof G. ; Grubmüller H. Constant pH Molecular Dynamics in Explicit Solvent with λ-Dynamics. J. Chem. Theory Comput. 2011, 7 , 1962–1978.21687785
(10) Goh G. B. ; Knight J. L. ; Brooks C. L. Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent. J. Chem. Theory Comput. 2012, 8 , 36–46.22337595
(11) Wallace J. A. ; Shen J. K. Charge-Leveling and Proper Treatment of Long-Range Electrostatics in All-Atom Molecular Dynamics at Constant pH. J. Chem. Phys. 2012, 137 , 184105.23163362
(12) Martins de Oliveira V. ; Liu R. ; Shen J. Constant pH Molecular Dynamics Simulations: Current Status and Recent Applications. Curr. Opin. Struct. Biol. 2022, 77 , 102498.36410222
(13) Huang Y. ; Harris R. C. ; Shen J. Generalized Born Based Continuous Constant pH Molecular Dynamics in Amber: Implementation, Benchmarking and Analysis. J. Chem. Inf. Model. 2018, 58 , 1372–1383.29949356
(14) Baptista A. M. ; Teixeira V. H. ; Soares C. M. Constant-pH Molecular Dynamics Using Stochastic Titration. J. Chem. Phys. 2002, 117 , 4184–4200.
(15) Goh G. B. ; Hulbert B. S. ; Zhou H. ; Brooks C. L. III Constant pH Molecular Dynamics of Proteins in Explicit Solvent with Proton Tautomerism: Explicit Solvent CPHMD of Proteins. Proteins 2014, 82 , 1319–1331.24375620
(16) Huang Y. ; Chen W. ; Wallace J. A. ; Shen J. All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water. J. Chem. Theory Comput. 2016, 12 , 5411–5421.27709966
(17) Radak B. K. ; Chipot C. ; Suh D. ; Jo S. ; Jiang W. ; Phillips J. C. ; Schulten K. ; Roux B. Constant-pH Molecular Dynamics Simulations for Large Biomolecular Systems. J. Chem. Theory Comput. 2017, 13 , 5933–5944.29111720
(18) Aho N. ; Buslaev P. ; Jansen A. ; Bauer P. ; Groenhof G. ; Hess B. Scalable Constant pH Molecular Dynamics in GROMACS. J. Chem. Theory Comput. 2022, 18 , 6148–6160.36128977
(19) Goh G. B. ; Knight J. L. ; Brooks III C. L. Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent. J. Chem. Theory Comput. 2012, 8 , 36–46.22337595
(20) Panahi A. ; Brooks C. L. Membrane Environment Modulates the p K a Values of Transmembrane Helices. J. Phys. Chem. B 2015, 119 , 4601–4607.25734901
(21) Shirts M. R. ; Pande V. S. Solvation Free Energies of Amino Acid Side Chain Analogs for Common Molecular Mechanics Water Models. J. Chem. Phys. 2005, 122 , 134508.15847482
(22) Hess B. ; van der Vegt N. F. A. Hydration Thermodynamic Properties of Amino Acid Analogues: A Systematic Comparison of Biomolecular Force Fields and Water Models. J. Phys. Chem. B 2006, 110 , 17616–17626.16942107
(23) Debiec K. T. ; Gronenborn A. M. ; Chong L. T. Evaluating the Strength of Salt Bridges: A Comparison of Current Biomolecular Force Fields. J. Phys. Chem. B 2014, 118 , 6561–6569.24702709
(24) Zhang H. ; Jiang Y. ; Cui Z. ; Yin C. Force Field Benchmark of Amino Acids. 2. Partition Coefficients between Water and Organic Solvents. J. Chem. Inf. Model. 2018, 58 , 1669–1681.30047730
(25) Zhang H. ; Yin C. ; Jiang Y. ; van der Spoel D. Force Field Benchmark of Amino Acids: I. Hydration and Diffusion in Different Water Models. J. Chem. Inf. Model. 2018, 58 , 1037–1052.29648448
(26) Ahmed M. C. ; Papaleo E. ; Lindorff-Larsen K. How Well Do Force Fields Capture the Strength of Salt Bridges in Proteins? PeerJ 2018, 6 , e4967.29910983
(27) Shen J. K. Uncovering Specific Electrostatic Interactions in the Denatured States of Proteins. Biophys. J. 2010, 99 , 924–932.20682271
(28) Harris J. A. ; Liu R. ; Martins de Oliveira V. ; Vázquez-Montelongo E. A. ; Henderson J. A. ; Shen J. GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber. J. Chem. Theory Comput. 2022, 18 , 7510–7527.36377980
(29) MacKerell A. D. ; Bashford D. ; Bellott M. ; Dunbrack R. L. ; Evanseck J. D. ; Field M. J. ; Fischer S. ; Gao J. ; Guo H. ; Ha S. ; Joseph-McCarthy D. ; Kuchnir L. ; Kuczera K. ; Lau F. T. K. ; Mattos C. ; Michnick S. ; Ngo T. ; Nguyen D. T. ; Prodhom B. ; Reiher W. E. ; Roux B. ; Schlenkrich M. ; Smith J. C. ; Stote R. ; Straub J. ; Watanabe M. ; Wiórkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102 , 3586–3616.24889800
(30) Mackerell A. D. ; Feig M. ; Brooks C. L. Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25 , 1400–1415.15185334
(31) Best R. B. ; Zhu X. ; Shim J. ; Lopes P. E. M. ; Mittal J. ; Feig M. ; MacKerell A. D. Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ 1 and χ 2 Dihedral Angles. J. Chem. Theory Comput. 2012, 8 , 3257–3273.23341755
(32) Huang J. ; Rauscher S. ; Nawrocki G. ; Ran T. ; Feig M. ; de Groot B. L. ; Grubmüller H. ; MacKerell A. D. CHARMM36m: An Improved Force Field for Folded and Intrinsically Disordered Proteins. Nat. Methods 2017, 14 , 71–73.27819658
(33) Maier J. A. ; Martinez C. ; Kasavajhala K. ; Wickstrom L. ; Hauser K. E. ; Simmerling C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11 , 3696–3713.26574453
(34) Tian C. ; Kasavajhala K. ; Belfon K. A. A. ; Raguette L. ; Huang H. ; Migues A. N. ; Bickel J. ; Wang Y. ; Pincay J. ; Wu Q. ; Simmerling C. ff19SB: Amino-Acid-Specific Protein Backbone Parameters Trained against Quantum Mechanics Energy Surfaces in Solution. J. Chem. Theory Comput. 2020, 16 , 528–552.31714766
(35) Sequeira J. G. N. ; Rodrigues F. E. P. ; Silva T. G. D. ; Reis P. B. P. S. ; Machuqueiro M. Extending the Stochastic Titration CpHMD to CHARMM36m. J. Phys. Chem. B 2022, 126 , 7870–7882.36190807
(36) Jorgensen W. L. ; Chandrasekhar J. ; Madura J. D. ; Impey R. W. ; Klein M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 , 926–935.
(37) Izadi S. ; Anandakrishnan R. ; Onufriev A. V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014, 5 , 3863–3871.25400877
(38) Alexov E. ; Mehler E. L. ; Baker N. ; Baptista M. , A.; Huang Y. ; Milletti F. ; Erik Nielsen J. ; Farrell D. ; Carstensen T. ; Olsson M. H. M. ; Shen J. K. ; Warwicker J. ; Williams S. ; Word J. M. Progress in the Prediction of p K a Values in Proteins: Prediction of p K a Values in Proteins. Proteins 2011, 79 , 3260–3275.22002859
(39) Hwang W. ; Austin S. L. ; Blondel A. ; Boittier E. D. ; Boresch S. ; Buck M. ; Buckner J. ; Caflisch A. ; Chang H.-T. ; Cheng X. ; Choi Y. K. ; Chu J.-W. ; Crowley M. ; Cui Q. ; Deng Y. ; Devereux M. ; Ding X. ; Feig M. ; Gao J. ; Glowacki D. R. ; Ii G. ; Hamaneh M. B. ; Harder E. D. ; Hayes R. L. ; Huang J. ; Huang Y. ; Im W. ; Islam S. M. ; Jiang W. ; Jones M. R. ; Kaser S. ; Kearns F. L. ; Kern N. R. ; Klauda J. B. ; Lazaridis T. ; Lee J. ; Lemkul J. A. ; Liu X. ; Luo Y. ; MacKerell A. D. ; Major D. T. ; Meuwly M. ; Nam K. ; Nilsson L. ; Ovchinnikov V. ; Paci E. ; Park S. ; Pastor R. W. ; Post C. B. ; Prasad S. ; Pu J. ; Qi Y. ; Rathinavelan T. ; Roe D. R. ; Roux B. ; Rowley C. N. ; Shen J. ; Simmonett A. C. ; Sodt A. J. ; Topfer K. ; Upadhyay M. CHARMM at 45: Comprehensive, Fast, and Free, a 15-Year Update. J. Phys. Chem. 2024, xxxx.
(40) Case D. A. ; Ben-Shalom I. Y. ; Brozell S. R. ; Cerutti D. S. ; Cheatham T. III ; Cruzeiro V. W. D. ; Darden T. A. ; Duke R. E. ; Ghoreishi D. ; Gilson M. K. ; Gohlke H. ; Goetz A. W. ; Greene D. ; Harris R. ; Homeyer N. ; Huang Y. ; Izadi S. ; Kovalenko A. ; Kurtzman T. ; Lee T. S. ; LeGrand S. ; Li P. ; Lin C. ; Liu J. ; Luchko T. ; Luo R. ; Mermelstein D. J. ; Merz K. M. ; Miao Y. ; Monard G. ; Nguyen C. ; Nguyen H. ; Omelyan I. ; Onufriev A. ; Pan F. ; Qi R. ; Roe D. R. ; Roitberg A. ; Sagui C. ; Schott-Verdugo S. ; Shen J. ; Simmerling C. L. ; Smith J. ; Salomon-Ferrer R. ; Swails J. ; Walker R. C. ; Wang J. ; Wei H. ; Wolf R. M. ; Wu X. ; Xiao L. ; York D. M. ; Kollman P. A. AMBER 2024. 2024.
(41) Luo Y. ; Roux B. Simulation of Osmotic Pressure in Concentrated Aqueous Salt Solutions. J. Phys. Chem. Lett. 2010, 1 , 183–189.
(42) Venable R. M. ; Luo Y. ; Gawrisch K. ; Roux B. ; Pastor R. W. Simulations of Anionic Lipid Membranes: Development of Interaction-Specific Ion Parameters and Validation Using NMR Data. J. Phys. Chem. B 2013, 117 , 10183–10192.23924441
(43) Yoo J. ; Aksimentiev A. Improved Parametrization of Li +, Na +, K +, and Mg 2+ Ions for All-Atom Molecular Dynamics Simulations of Nucleic Acid Systems. J. Phys. Chem. Lett. 2012, 3 , 45–50.
(44) Yoo J. ; Aksimentiev A. Improved Parameterization of Amine–Carboxylate and Amine–Phosphate Interactions for Molecular Dynamics Simulations Using the CHARMM and AMBER Force Fields. J. Chem. Theory Comput. 2016, 12 , 430–443.26632962
(45) Yoo J. ; Aksimentiev A. New Tricks for Old Dogs: Improving the Accuracy of Biomolecular Force Fields by Pair-Specific Corrections to Non-Bonded Interactions. Phys. Chem. Chem. Phys. 2018, 20 , 8432–8449.29547221
(46) Ferguson N. ; Sharpe T. D. ; Schartau P. J. ; Sato S. ; Allen M. D. ; Johnson C. M. ; Rutherford T. J. ; Fersht A. R. Ultra-Fast Barrier-limited Folding in the Peripheral Subunit-binding Domain Family. J. Mol. Biol. 2005, 353 , 427–446.16168437
(47) Arbely E. ; Rutherford T. J. ; Sharpe T. D. ; Ferguson N. ; Fersht A. R. Downhill versus Barrier-Limited Folding of BBL 1: Energetic and Structural Perturbation Effects upon Protonation of a Histidine of Unusually Low pKa. J. Mol. Biol. 2009, 387 , 986–992.19136007
(48) Arbely E. ; Rutherford T. J. ; Neuweiler H. ; Sharpe T. D. ; Ferguson N. ; Fersht A. R. Carboxyl pKa Values and Acid Denaturation of BBL. J. Mol. Biol. 2010, 403 , 313–327.20816989
(49) Joung I. S. ; Cheatham T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112 , 9020–9041.18593145
(50) Jorgensen W. L. ; Chandrasekhar J. ; Madura J. D. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 , 926.
(51) Beglov D. ; Roux B. Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations. J. Chem. Phys. 1994, 100 , 9050–9063.
(52) Lev B. ; Roux B. ; Noskov S. Y. Relative Free Energies for Hydration of Monovalent Ions from QM and QM/MM Simulations. J. Chem. Theory Comput. 2013, 9 , 4165–4175.26592407
(53) Nguyen H. ; Roe D. R. ; Simmerling C. Improved Generalized Born Solvent Model Parameters for Protein Simulations. J. Chem. Theory Comput. 2013, 9 , 2020–2034.25788871
(54) Piana S. ; Robustelli P. ; Tan D. ; Chen S. ; Shaw D. E. Development of a Force Field for the Simulation of Single-Chain Proteins and Protein–Protein Complexes. J. Chem. Theory Comput. 2020, 16 , 2494–2507.31914313
(55) Izadi S. ; Onufriev A. V. Accuracy Limit of Rigid 3-Point Water Models. J. Chem. Phys. 2016, 145 , 074501.27544113
(56) Zhang H. ; Jiang Y. ; Yan H. ; Yin C. ; Tan T. ; Van Der Spoel D. Free-Energy Calculations of Ionic Hydration Consistent with the Experimental Hydration Free Energy of the Proton. J. Phys. Chem. Lett. 2017, 8 , 2705–2712.28561580
(57) Nguyen H. ; Maier J. ; Huang H. ; Perrone V. ; Simmerling C. Folding Simulations for Proteins with Diverse Topologies Are Accessible in Days with a Physics-Based Force Field and Implicit Solvent. J. Am. Chem. Soc. 2014, 136 , 13959–13962.25255057
(58) Henderson J. A. ; Lui R. ; Harris J. A. ; Huang Y. ; de Oliviera V. M. ; Shen J. A Guide to the Continuous Constant pH Molecular Dynamics Methods in Amber and CHARMM [Article v1.0]. Liv. J. Comput. Mol. Sci. 2022, 4 , 1563.
(59) Henderson J. A. ; Verma N. ; Harris R. C. ; Liu R. ; Shen J. Assessment of Proton-Coupled Conformational Dynamics of SARS and MERS Coronavirus Papain-like Proteases: Implication for Designing Broad-Spectrum Antiviral Inhibitors. J. Chem. Phys. 2020, 153 , 115101.32962355
(60) Thurlkill R. L. ; Grimsley G. R. ; Scholtz J. M. ; Pace C. N. pK Values of the Ionizable Groups of Proteins. Protein Sci. 2006, 15 , 1214–1218.16597822
