Download Hi-Res ImageDownload to MS-PowerPointCite This:J. Chem. Theory Comput. 2022, 18, 10, 6148-6160

Molecular MechanicsSeptember 21, 2022

Scalable Constant pH Molecular Dynamics in GROMACS
Click to copy article linkArticle link copied!

Noora Aho*
Noora Aho
Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014Jyväskylä, Finland
*Email: [email protected]
More by Noora Aho
Pavel Buslaev*
Pavel Buslaev
Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014Jyväskylä, Finland
*Email: [email protected]
More by Pavel Buslaev
https://orcid.org/0000-0003-2031-4691
Anton Jansen
Anton Jansen
Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44Stockholm, Sweden
More by Anton Jansen
Paul Bauer
Paul Bauer
Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44Stockholm, Sweden
More by Paul Bauer
https://orcid.org/0000-0002-2268-0065
Gerrit Groenhof*
Gerrit Groenhof
Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014Jyväskylä, Finland
*Email: [email protected]
More by Gerrit Groenhof
https://orcid.org/0000-0001-8148-5334
Berk Hess*
Berk Hess
Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44Stockholm, Sweden
*Email: [email protected]
More by Berk Hess

Open PDF Supporting Information (3)

Journal of Chemical Theory and Computation

Cite this: J. Chem. Theory Comput. 2022, 18, 10, 6148–6160

Click to copy citationCitation copied!

https://doi.org/10.1021/acs.jctc.2c00516

Published September 21, 2022

CC-BY 4.0 .

Abstract

Click to copy section linkSection link copied!

Molecular dynamics (MD) computer simulations are used routinely to compute atomistic trajectories of complex systems. Systems are simulated in various ensembles, depending on the experimental conditions one aims to mimic. While constant energy, temperature, volume, and pressure are rather straightforward to model, pH, which is an equally important parameter in experiments, is more difficult to account for in simulations. Although a constant pH algorithm based on the λ-dynamics approach by Brooks and co-workers [Kong, X.; Brooks III, C. L. J. Chem. Phys.1996, 105, 2414–2423] was implemented in a fork of the GROMACS molecular dynamics program, uptake has been rather limited, presumably due to the poor scaling of that code with respect to the number of titratable sites. To overcome this limitation, we implemented an alternative scheme for interpolating the Hamiltonians of the protonation states that makes the constant pH molecular dynamics simulations almost as fast as a normal MD simulation with GROMACS. In addition, we implemented a simpler scheme, called multisite representation, for modeling side chains with multiple titratable sites, such as imidazole rings. This scheme, which is based on constraining the sum of the λ-coordinates, not only reduces the complexity associated with parametrizing the intramolecular interactions between the sites but also is easily extendable to other molecules with multiple titratable sites. With the combination of a more efficient interpolation scheme and multisite representation of titratable groups, we anticipate a rapid uptake of constant pH molecular dynamics simulations within the GROMACS user community.

This publication is licensed under

CC-BY 4.0 .

License Summary*
You are free to share(copy and redistribute) this article in any medium or format and to adapt(remix, transform, and build upon) the material for any purpose, even commercially within the parameters below:
- Creative Commons (CC): This is a Creative Commons license.
- Attribution (BY): Credit must be given to the creator.
View full license
*Disclaimer
This summary highlights only some of the key features and terms of the actual license. It is not a license and has no legal value. Carefully review the actual license before using these materials.
License Summary*
You are free to share(copy and redistribute) this article in any medium or format and to adapt(remix, transform, and build upon) the material for any purpose, even commercially within the parameters below:
- Creative Commons (CC): This is a Creative Commons license.
- Attribution (BY): Credit must be given to the creator.
View full license
*Disclaimer
This summary highlights only some of the key features and terms of the actual license. It is not a license and has no legal value. Carefully review the actual license before using these materials.

Subjects

Introduction

Click to copy section linkSection link copied!

Since their introduction more than four decades ago, molecular dynamics (MD) computer simulations have come of age. (1) Thanks to improvements in computer hardware, algorithmic developments, as well as increased accuracy of force fields, MD simulation has evolved into a predictive technique that can complement experiments by providing atomistic insights into the dynamics of complex systems. (1,2) While many experimental conditions can be modeled with good accuracy, the aqueous proton concentration, or pH, is typically accounted for indirectly by constraining the protonation states of titratable residues to their, presumed, most probable form at the start of the simulation. Because the electrostatic interactions depend critically on the protonation state of the residues, the pH affects the conformational ensemble. Conversely, because the conformation can influence the proton affinity of the residues, or pK_a, a direct correlation exists between pH and conformational dynamics, which cannot be captured if the protonation state is kept fixed in the simulation. (3)

To overcome this limitation in classical MD simulations and include the effect of pH on the conformational sampling directly, several solutions have been proposed in the last decades (4,5) and used to investigate pH-dependent protein–protein (6) and protein-RNA interactions, (7) drug binding, (8,9) and structural changes. (10,11) These solutions can be roughly divided into a category that relies on discrete changes in protonation states (12−19) and a second category in which a protonation state can change continuously. (20−33) More recently, a third category that relies on the transfer of proton-like particles between titratable sites, including protein residues and solvent molecules, was proposed for the Martini force field. (34)

In the discrete constant pH approaches, the protonation state of a residue can change at regular intervals of the simulation according to a Metropolis Monte Carlo criterion. (16,17,28,35,36) To avoid a low acceptance rate due to unfavorable solvent configurations, the Monte Carlo step is performed based on free energies calculated using the approximation of either an implicit solvent representation, (13,15) or a short all-atom thermodynamic integration. (14,37)

Most continuous approaches for MD at constant pH are based on the λ-dynamics technique developed by Brooks and co-workers. (38) A one-dimensional λ-coordinate with fictitious mass m_λ is introduced for each titratable site, and the equations of motion for these additional degrees of freedom are integrated along with the Cartesian positions of the atoms. (21) The λ-coordinate defines the protonation state of the residue: at λ = 0 the residue is protonated and interacts with the rest of the system as such, while at λ = 1 it is deprotonated. The energy function that acts on the λ-coordinates depends on (i) the intrinsic proton affinity (reference pK_a) of the titratable site in water, (ii) the interactions with the environment, which are mostly electrostatic, (39) and (iii) the pH of the solvent, which is set by the user. In addition, potentials are introduced to bias sampling toward the physical states at λ = 0 and λ = 1. Protons are not transferred directly between the titratable residues and the solvent molecules but rather exchanged with an external proton bath. Because the chemical potential of this bath is determined by the proton concentration (pH) of the aqueous solution, constant pH MD (CpHMD) simulations based on λ-dynamics are performed in a grand canonical ensemble for the proton degrees of freedom.

While λ-dynamics-based constant pH approaches were originally developed for implicit solvent simulations, (21) they have since then been adapted for explicit solvent simulations as well. (18,23−26,28,30) The key computational challenge for explicit solvent implementations is the long-range electrostatic interaction, for which multiple solutions have been suggested, including a shifted cutoff scheme, (26) a hybrid scheme combining the particle mesh Ewald (PME) treatment for the Cartesian coordinates with the generalized Born model for the λ-particles, (24,40) and a fully consistent PME treatment for both λ and Cartesian degrees of freedom. (23,30)

In addition to accurate modeling of the long-range electrostatic interactions, also sampling can pose a serious challenge to simulations at constant pH. While the choice for the PME method in the original implementation of λ-dynamics in the fork of GROMACS 3.3 release (23) was motivated by its accurate description of long-range electrostatics, the linear increase of the computational effort with the number of titratable sites limited the sampling efficiency, which meant that systems with many titratable sites could not be studied in practice.

To remove this bottleneck and enable constant pH MD with GROMACS at a modest additional cost compared to a standard simulation, we switch to an alternative scheme for computing the long-range electrostatic interactions of the λ-particles under periodic boundary conditions. The alternative scheme is based on a linear interpolation of partial charges (21) rather than the potential energy functions as in the original implementation of constant pH MD in a GROMACS fork. (23)

Although the previous implementation of constant pH in a GROMACS fork was documented and shared with the community as an open-source program, there has been some misunderstanding about how electrostatic interactions were computed for λ-particles. (24,25,30) To resolve this controversy, we first explain in detail how the electrostatic interactions were calculated in the previous GROMACS implementation of constant pH MD. We next contrast this linear interpolation between the potential energy functions of the protonated and deprotonated states of a residue on the one hand to the interpolation between the partial charges of both states on the other hand (21) and show why the latter is computationally much more efficient. We then demonstrate the superior performance of the charge-interpolation scheme by running a series of constant pH MD simulations of amino acids and proteins. To emphasize that the new constant pH implementation in GROMACS is not restricted to a specific force field (or the resolution of a force field model) nor to a specific algorithm for evaluating electrostatic interactions, we also show the results of constant pH MD simulations with the Martini coarse-grained force field, (41) in combination with a shifted cutoff electrostatics model. Because of GROMACS’ large user community, we expect our work to increase the popularity of constant pH MD simulations.

Theory

Click to copy section linkSection link copied!

Before discussing the differences between linear interpolating of the potential energy functions on the one hand, (23) and of partial charges on the other hand, (21) for computing the potential energy landscape of the titration coordinates, we briefly review the λ-dynamics approach (38) that forms the basis for the constant pH molecular dynamics algorithm in GROMACS. (23)

λ-Dynamics-Based Constant pH MD Simulations

A titratable site i can exist in a protonated or deprotonated state. The protonation state affects the interactions between the site and the rest of the system. In constant pH MD simulations based on λ-dynamics, (38) an additional coordinate λ_i is introduced for each site i, and the potential energy function of the total system is continuously interpolated between the two protonation states along this coordinate, i.e., V(λ_i) . (21) A fictitious mass m_λ is assigned to each λ_i-coordinate, and the coordinates evolve along with the Cartesian degrees of freedom of all atoms in the system, based on Newton’s equations of motion. Thus, the total Hamiltonian of the system is

H (R, λ) = \sum_{i}^{N_{sites}} \frac{m_{λ}}{2} {\dot{λ}}_{i}^{2} + \sum_{j}^{N_{atoms}} \frac{m_{j}}{2} {\dot{r}}_{j}^{2} + V (R, λ)

(1)

where R is the vector of the Cartesian coordinates r_j of all N_atoms atoms with mass m_j and λ is the vector of the λ_i coordinates of all N_sites titratable sites.

λ-Dependent Potential Energy Function

In addition to the interpolation between the potentials of the protonated V_A(R) and deprotonated states V_B(R), three more λ-dependent terms are included in the potential energy function of the total system V(R, λ), as illustrated in Figure 1: (i) a correction potential V_i^MM(λ_i) to compensate for missing quantum mechanical contributions to proton affinities (Figure 1B); (ii) a biasing potential V_i^bias(λ_i) that enhances sampling of the physical end states at λ_i = 0 and λ_i = 1 (Figure 1C); and (iii) a pH-dependent term V^pH(λ_i) to model the chemical potential of protons in water (Figure 1D).

Figure 1. Illustration of the additional λ-dependent potential energy terms (B–D). Panel A shows the protonation free energy of a titratable residue in its reference state obtained at the force field level, ΔG_i^MM(λ_i) . To compensate for the shortcomings of the force field and obtain a zero free energy difference between the two protonation states (A + B), we add the correction potential, V_i^MM(λ), shown in panel B. A biasing potential, V_i^bias(λ), (42) is introduced to avoid sampling of nonphysical states (panel C). To model the proton chemical potential (pH), we add a pH-dependent term, V^pH(λ_i) (panel D). For a titratable residue at pH ≠ pK_a, the total λ-dependent potential, including the interpolated force field functions and the three additional terms, is illustrated in panel (A + B + C + D).

The purpose of adding the correction term V_i^MM(λ_i) (Figure 1B) is to make the interpolated potential function flat if the titratable site i is in its reference state, for which the proton affinity is known experimentally, at pH = pK_a,i. This potential is determined by evaluating the deprotonation free energy of the single residue in water (reference state) at the force field level by thermodynamic integration along the λ-coordinate (Figure 1A):

V_{i}^{MM} (λ_{i}) = - Δ G_{i}^{MM} (λ_{i})

(2)

To prevent sampling of the nonphysical states between λ_i = 0 and λ_i = 1 on this flat potential energy surface while still enabling sufficient transitions between the physical end-states to sample both protonation states with the correct thermodynamic weight, we introduce the biasing potential V_i^bias(λ_i) suggested by Donnini et al. (42) (Figure 1C).

The pH-dependent term V^pH(λ_i) (Figure 1D) is a correction that includes the effect of the solution pH on the free energy difference between the protonated and deprotonated states, such that this difference is

V^{pH} (λ_{i} = 1) - V^{pH} (λ_{i} = 0) = R T \ln (10) [p K_{a, i} - pH]

(3)

where we use the experimentally determined pK_a,i value of residue i in its reference state. Although various forms for this potential have been suggested, (29,42,43) we propose a smooth step-function-based potential:

\begin{matrix} V^{pH} (λ_{i}) & = & R T \ln (10) [p K_{a, i} - pH] \times \frac{1}{1 + \exp (- 2 k_{1} (λ_{i} - 1 + x_{0}))}, if pH > p K_{a} \\ V^{pH} (λ_{i}) & = & R T \ln (10) [p K_{a, i} - pH] \times \frac{1}{1 + \exp (- 2 k_{1} (λ_{i} - x_{0}))}, if pH \leq p K_{a} \end{matrix}

(4)

where k₁ and x₀ define the steepness and kink position of the step function. In this form, illustrated in Figure 1D, the pH-dependent potential also aids in preventing the sampling of nonphysical states, i.e., 0.1 < λ < 0.9.

Linear Interpolation of Potential Energy Functions

In the previous implementation of constant pH MD in a GROMACS fork, the smooth interpolation of the force field potential energy function between the protonated and deprotonated states was achieved by linearly interpolating the force field potentials of these states. (23) Thus, for a single titratable site, the λ-dependent potential is given by

V (R, λ) = (1 - λ) V_{A} (R) + λ V_{B} (R) + V^{MM} (λ) + V^{bias} (λ) + V^{pH} (λ)

(5)

with V^MM(R, λ), V^bias(λ), and V^pH(λ) the correction, biasing, and pH-dependent potentials, respectively, that were briefly discussed above, and with short-hand notations for

\begin{matrix} V_{A} (R) = V (R, λ = 0) \\ V_{B} (R) = V (R, λ = 1) \end{matrix}

The gradient required for updating λ according to Newton’s equations of motion is

\frac{\partial V (R, λ)}{\partial λ} = V_{B} (R) - V_{A} (R) + \frac{\partial}{\partial λ} V^{MM} (λ) + \frac{\partial}{\partial λ} V^{bias} (λ) + \frac{\partial}{\partial λ} V^{pH} (λ)

(6)

Thus, the evaluation of the force on the λ-particle requires that the potential energy, including the long-range electrostatic interactions, is computed twice: once for λ = 0 (i.e., V_A) and once more for λ = 1 (i.e., V_B). If the Particle-Mesh-Ewald (PME) method is used to compute those long-range electrostatic interactions, (44,45) separate PME grid builds are needed because the charge distributions are not identical in states A and B.

For systems with many titratable sites, multiple λ-groups are introduced. Because the analytical expressions for the correction, biasing, and pH-dependent terms in eq 5 are additive, we no longer consider them explicitly in what follows and focus exclusively on the interpolation of the force field potential energies between the multiple protonation states of the system. For N λ-coordinates, there are 2^N such states and the interpolation generalizes to (20)

V (R, λ) = \sum_{l}^{2^{N}} V (R, l) Π_{i}^{N} λ_{i}^{(l_{i})} {(1 - λ_{i})}^{1 - l_{i}}

(7)

Here, we represent the N λ_i-coordinates as an N-dimensional vector λ. The 2^N possible protonation states of the system are represented by the N-dimensional vector l with elements l_i equal to 0 or 1 that indicate whether a site i is protonated (λ_i = 0) or deprotonated (λ_i = 1). The sum runs over all 2^N possible combinations of l_i = 0 and l_i = 1. The gradient required for updating λ_i is obtained by deriving the interpolated potential, V(R, λ), with respect to λ_i:

(8)

where the ′ indicates that λ_i is omitted from vector λ. Note that, as we focus only on the interpolated potentials, the biasing, correction, and pH-dependent terms are left out.

In general, the number of terms in the potential (eq 7) increases exponentially with the number of titratable sites. However, for pairwise interactions involving titratable sites whose nonbonded force field parameters do not depend on the protonation state of the other sites (chemically uncoupled sites), the number of terms required to evaluate the interpolated potential scales linearly. For systems with such “chemically”, or “topologically” uncoupled sites, the interpolated potential contains four types of interactions

V (R, λ) = V^{rest - rest} (R) + V^{λ - rest} (R, λ) + V^{λ - λ} (R, λ) + V^{λ} (R, λ)

(9)

For pairwise electrostatic interactions, the terms on the right-hand side are defined as

1.	Interactions between atoms that are not part of any λ-group, and hence independent of the λ_i’s: $V^{rest - rest} (R) = \frac{1}{2} \sum_{i}^{n_{rest}} \sum_{j}^{n_{rest}} \frac{q_{i} q_{j}}{4 π ϵ_{0} r_{i j}}$ (10) where the sums run over all n_rest atoms that are not part of a λ-group.
2.	Interpolated interactions between atoms of each λ-group with atoms that are not part of any λ-group: $V^{λ - rest} (R, λ) = \sum_{k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{rest}} {(1 - λ_{k}) \frac{q_{i}^{A} q_{j}}{4 π ϵ_{0} r_{i j}} + λ_{k} \frac{q_{i}^{B} q_{j}}{4 π ϵ_{0} r_{i j}}}$ (11) where the first sum runs over all titratable sites, the second one runs over all n_k atoms of the k th λ-group, and the final sum runs over all atoms that are not part of any λ group.
3.	Interpolated interactions between atoms belonging to two different λ-groups: $V^{λ - λ} (R, λ) = \sum_{k}^{N_{sites}} {(1 - λ_{k}) \sum_{m, m \neq k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} [(1 - λ_{m}) \frac{q_{i}^{A} q_{j}^{A}}{4 π ϵ_{0} r_{i j}} + λ_{m} \frac{q_{i}^{A} q_{j}^{B}}{4 π ϵ_{0} r_{i j}}] + λ_{k} \sum_{m, m \neq k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} [(1 - λ_{m}) \frac{q_{i}^{B} q_{j}^{A}}{4 π ϵ_{0} r_{i j}} + λ_{m} \frac{q_{i}^{B} q_{j}^{B}}{4 π ϵ_{0} r_{i j}}]} = \sum_{k}^{N_{sites}} \sum_{m, m \neq k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} \frac{[(1 - λ_{k}) q_{i}^{A} + λ_{k} q_{i}^{B}] [(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}}$ (12)
4.	Interpolated interactions within each of the λ-groups: $V^{λ} (R, λ) = \frac{1}{2} \sum_{k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{k}} {(1 - λ_{k}) \frac{q_{i}^{A} q_{j}^{A}}{4 π ϵ_{0} r_{i j}} + λ_{k} \frac{q_{i}^{B} q_{j}^{B}}{4 π ϵ_{0} r_{i j}}}$ (13)

From eq 8, the gradient with respect to λ_k is

\begin{matrix} \frac{\partial}{\partial λ_{k}} V_{Coul} (R, λ) & = & V_{Coul} (R, λ^{'}, λ_{k} = 1) - V_{Coul} (R, λ^{'}, λ_{k} = 0) \\ = & \sum_{i}^{n_{k}} \sum_{j}^{n_{rest}} [\frac{q_{i}^{B} q_{j}}{4 π ϵ_{0} r_{i j}} - \frac{q_{i}^{A} q_{j}}{4 π ϵ_{0} r_{i j}}] + \sum_{m, m \neq k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} [\frac{q_{i}^{B} [(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}} - \frac{q_{i}^{A} [(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}}] + \frac{1}{2} \sum_{i}^{n_{k}} \sum_{j}^{n_{k}} [\frac{q_{i}^{B} q_{j}^{B}}{4 π ϵ_{0} r_{i j}} - \frac{q_{i}^{A} q_{j}^{A}}{4 π ϵ_{0} r_{i j}}] \end{matrix}

(14)

Thus, the evaluation of the Coulomb contribution to the gradient for each λ_k-group requires two electrostatic computations, with the interpolated partial charges of the other λ_m sites (i.e., q_j(λ_m) = (1 – λ_m)q_j^A + λ_mq_j^B):

\begin{matrix} \frac{\partial}{\partial λ_{k}} V_{Coul} (R, λ) & = & \sum_{i}^{n_{k}} q_{i}^{B} \times [\sum_{a}^{n_{rest}} \frac{q_{a}}{4 π ϵ_{0} r_{i a}} + \sum_{m, m \neq k}^{N_{sites}} \sum_{j}^{n_{m}} \frac{(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}}{4 π ϵ_{0} r_{i j}} + \frac{1}{2} \sum_{j}^{n_{k}} \frac{q_{j}^{B}}{4 π ϵ_{0} r_{i j}}] - \sum_{i}^{n_{k}} q_{i}^{A} \times [\sum_{a}^{n_{rest}} \frac{q_{a}}{4 π ϵ_{0} r_{i a}} + \sum_{m, m \neq k}^{N_{sites}} \sum_{j}^{n_{m}} \frac{(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}}{4 π ϵ_{0} r_{i j}} + \frac{1}{2} \sum_{j}^{n_{k}} \frac{q_{j}^{A}}{4 π ϵ_{0} r_{i j}}] \\ = & \sum_{i}^{n_{k}} q_{i}^{B} Φ_{k}^{B} (R_{i}, λ^{'}) - \sum_{i}^{n_{k}} q_{i}^{A} Φ_{k}^{A} (R_{i}, λ^{'}) \end{matrix}

(15)

Here, we introduced the electrostatic potential Φ_k^A(R_i, λ′) of the system with partial charges of λ-group k in the protonated state (q_i^A) and interpolated charges for all other λ-groups. As before, λ′ is the vector with all λ_m’s except λ_k. Likewise, electrostatic potential Φ_k^B(R_i, λ′) is evaluated with the partial charges of λ-group k in the deprotonated state (q_i^B) and the same interpolated charges for all other λ-groups. Thus, 2N_sites computations are needed to evaluate the gradients for all titratable sites. The same arguments apply to pairwise Lennard-Jones interactions, but because the contribution of Lennard-Jones interaction to pK_a shift is minor, we neglected them in this work (see Supporting Information).

Linear Interpolation of Partial Charges

While the linear scaling of the gradients for the pairwise potentials in eq 15 in principle is a great improvement over the formal exponential scaling in eq 8, the requirement of performing 2N_sites calculations per MD step still poses a computational bottleneck, in particular for larger systems. To overcome this bottleneck for electrostatic interactions, we follow the suggestion by Brooks and co-workers to interpolate charges rather than interaction functions. (21) When interpolating the partial charges between the protonation states of N_sites chemically uncoupled titratable sites, the λ-dependent Coulomb energy becomes

V_{coul} (R, λ) = V^{rest - rest} (R) + V^{λ - rest} (R, λ) + V^{λ - λ} (R, λ) + V^{λ} (R, λ) = \frac{1}{2} \sum_{i}^{n_{rest}} \sum_{j}^{n_{rest}} \frac{q_{i} q_{j}}{4 π ϵ_{0} r_{i j}} + \sum_{k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{i}^{n_{rest}} \frac{((1 - λ_{k}) q_{i}^{A} + λ_{k} q_{i}^{B}) q_{j}}{4 π ϵ_{0} r_{i j}} + \sum_{k}^{N_{sites}} \sum_{m; m \neq k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} \frac{[(1 - λ_{k}) q_{i}^{A} + λ_{k} q_{i}^{B}] [(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}} + \frac{1}{2} \sum_{k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{k}} \frac{[(1 - λ_{k}) q_{i}^{A} + λ_{k} q_{i}^{B}] [(1 - λ_{k}) q_{j}^{A} + λ_{k} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}}

(16)

The gradient of the potential energy with respect to λ_k is

\frac{\partial V_{coul} (R, λ)}{\partial λ_{k}} = \sum_{i}^{n_{k}} \sum_{j}^{n_{rest}} \frac{(q_{i}^{B} - q_{i}^{A}) q_{j}}{4 π ϵ_{0} r_{i j}} + \sum_{m = 1; m \neq k}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} \frac{(q_{i}^{B} - q_{i}^{A}) [(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}} + \sum_{i}^{n_{k}} \sum_{j}^{n_{k}} \frac{(q_{i}^{B} - q_{i}^{A}) [(1 - λ_{k}) q_{j}^{A} + λ_{k} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}} = \sum_{i}^{n_{k}} \sum_{j}^{n_{rest}} \frac{(q_{i}^{B} - q_{i}^{A}) q_{j}}{4 π ϵ_{0} r_{i j}} + \sum_{m = 1}^{N_{sites}} \sum_{i}^{n_{k}} \sum_{j}^{n_{m}} \frac{(q_{i}^{B} - q_{i}^{A}) [(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}]}{4 π ϵ_{0} r_{i j}} = \sum_{i}^{n_{k}} Φ (R_{i}, λ) Δ q_{i}

(17)

where Φ(R_i, λ) is the electrostatic potential at the position of atom i due to the charge distribution of all other atoms in the system, including the atoms of all titratable sites, for which the partial charges are interpolated:

Φ (R_{i}, λ) = \sum_{j}^{n_{rest}} \frac{q_{j}}{4 π ϵ_{0} r_{i j}} + \sum_{m}^{N_{sites}} \sum_{j}^{n_{m}} \frac{(1 - λ_{m}) q_{j}^{A} + λ_{m} q_{j}^{B}}{4 π ϵ_{0} r_{i j}}

and Δq_i is the difference between the atomic charges of titratable residue i in the protonated (A) and deprotonated (B) states:

Δ q_{i} = q_{i}^{B} - q_{i}^{A}

In contrast to when potential energy functions are interpolated, the same electrostatic potential is used to evaluate the electrostatic forces on both the atoms and the λ-particles. Therefore, a single electrostatic calculation per time step suffices. If the electrostatic interactions are modeled with the smooth Particle Mesh Ewald method, (44,45) the short-range real-space interactions and long-range reciprocal-space interactions are computed separately. For the pairwise short-range interactions in real space, an additional calculation for each interacting pair and a subsequent accumulation of the potential at each atom is needed. Whereas this calculation comes at no extra computational cost if the standard pair interaction kernels are used, the accumulation leads to a measurable computational overhead, as we will show later. For the mesh part of the PME calculation, a gathering of potentials from the grid is required for charges in λ-groups only, but this also comes at a negligible computational overhead. Because the extra effort required to compute the gradients on the λ-particles is rather small, a constant pH MD implementation based on charge interpolation is computationally not much more expensive than a normal MD simulation, which is a major improvement with respect to the previous CpHMD implementation in GROMACS. (23)

Multisite Representation of Chemically Coupled Titratable Sites

If titratable sites are “chemically” or “topologically” coupled, the force field parameters of one site depend on the value of the λ-coordinate of the other site, and vice versa. For example, histidine can exist in three protonation states, as shown in Figure 2. In most force fields, the partial charges of all atoms in the His side chain, including the two sites, depend on the protonation state. Hence, if the N_δ site changes protonation, the electrostatic interactions of the N_ϵ site are also affected.

Figure 2. Multisite representation illustrated for a histidine side chain. Each possible protonation state is represented by its own λ-coordinate. HSP refers to doubly protonated histidine, HSD, and HSE to histidine singly protonated at the N_δ or N_ϵ, respectively. HS0 is a common, nonphysical state of the residue. To restrict the sampling to the plane connecting the physical states, a constraint λ₁ + λ₂ + λ₃ = 1 is applied (gray plane). A biasing potential is also applied to enhance sampling at the end states, where one of the λ-coordinates is close to one, while the other coordinates are close to zero.

To model the chemically coupled sites in the histidine side chain, Khandogin and Brooks introduced two λ-coordinates: (22) one that interpolates between the double and single protonated forms and a second coordinate switching between protonation at the N_δ and the N_ϵ atoms. Donnini et al. introduced separate λ-coordinates for N_δ and N_ϵ. (23) In both solutions, the coupling between the coordinates is achieved with a two-dimensional correction potential.

Because extending the dimensionality beyond two coordinates is difficult from both the implementation and parametrization perspective, Brooks and co-workers introduced a multisite representation, (46,47) where a separate λ_i,k-coordinate is assigned to each physical state k of a titratable group i. For a residue with multiple “chemically”-coupled titratable sites, each λ_i,k-coordinate has the same state at λ_i,k = 0, while at λ_i,k = 1, the group is in one of the n_i possible protonation states (i.e., state k) of residue i. The state at λ_i,k = 0 is the same for all λ_i,k-coordinates in residue i but does not correspond to a physical protonation state of the residue and neither do states for which the sum of the λ_i,k-s is not equal to 1 (Figure 2). To restrict sampling to the (hyper-)plane connecting the physical states, the sum of the λ_i,k values is constrained (∑_kλ_i,k = 1). Since the λ-dynamics implementation in a fork of GROMACS relies on linear λ-coordinates, rather than on auxiliary circular coordinates that would fulfill the constraints by construction, (46,48) we apply a constraint on the sum of λ_i,k-coordinates. To efficiently apply this constraint, we use an analytical expression to solve a generalized version of the charge constraint introduced by Donnini et al. (42) (see the Appendix). While an atom can be part of multiple λ_i,k-coordinates in residue i, each affecting its charge, we show in the Supporting Information that the expression for the contribution of this atom to the total Coulomb energy is identical to that of an uncoupled site (eq 17).

In the multisite representation, each λ-coordinate is independent of the others and thus evolves on a one-dimensional potential (eq 4), similar to that of “chemically” uncoupled sites. However, in contrast to the uncoupled sites, the correction potential V^MM is multidimensional as its value depends on all λ_i,k-coordinates representing each of the possible protonation states of residue i. These potentials are obtained through a least-squares fit of a multidimensional polynomial to the ensemble-averaged gradients of the potentials with respect to λ_i,k evaluated on the (n_i – 1)-dimensional grid of the n_i coupled λ-coordinates, i.e., ⟨∂ V/∂λ_i,k⟩_{λ₁...λ_ni}. The fitting procedure is explained in detail in the Supporting Information.

The multisite representation can be applied to residues with any number of titratable sites, including residues with only a single titratable site. In the latter case, two λ-coordinates, corresponding to the protonated state (λ_i,1 = 1, λ_i,2 = 0) and deprotonated state (λ_i,1 = 0, λ_i,2 = 1), are introduced with a constraint on their sum (λ_i,1 + λ_i,2 = 1).

Methods

Click to copy section linkSection link copied!

We have implemented the algorithms for CpHMD with charge interpolation in a fork of GROMACS software package (2021 release). (49) The code and manuals are available for free at https://gitlab.com/gromacs-constantph/constantph. Here we verify the validity of our implementation for reproducing pK_a values of peptides and proteins. To demonstrate that the linear interpolation of charges (eq 17) scales better with the number of titratable sites in the system than the linear interpolation of interaction functions (eq 15), we compared the scaling between our new implementation, which is based on linear charge interpolation on the one hand, and a previous implementation in a fork of the GROMACS 3.3 release, which is based on linear interpolation of the force field potentials on the other hand. (23) To estimate the additional computational effort required for performing CpHMD with the new implementation, we also compared the performance of a CpHMD simulation to that of a normal MD simulations on both CPUs and GPUs.

Simulated Systems

To test the implementation, we performed constant pH MD simulations of the following systems: (1) glutamic acid (Glu), (2) aspartic acid (Asp), (3) histidine (His), (4) Cardiotoxin V (PDB ID: 1CVO (50)), (5) hen egg white lyzozyme (HEWL, PDB ID: 2LZT (51)), (6) the GLIC pentameric ligand-gated ion channel (PDB ID: 4HFI (52)), and (7) turkey ovomucoid inhibitor (PDB ID: 2GKR (53)).

In systems 1–6, the interactions were modeled with the CHARMM36 (54) all-atom (AA) force field, with some modifications in the torsion parameters to accelerate the convergence. These modifications are presented and validated in an accompanying paper, in which we report the application of our constant pH implementation to lysine, C-, and N-termini. (55) Systems 1–5 were also simulated with the Martini 2.0 coarse grained (CG) force field. (41) The martinize.py script was used to automatically generate the CG representation of these systems. (56) System 7 was simulated to compare the efficiency of interpolating charges and potentials. The interactions in this system were modeled with the OPLS force field (57) because the GROMACS 3.3 release, on which the linear interpolation of potentials implementation was based, does not support the CMAP correction that is needed for the CHARMM36 force field. (58)

The amino acids Glu, Asp, and His were modeled as tripeptides Ala-X-Ala with acetylated N-terminus (ACE) and N-methylamidated C-terminus (CT₃). The proteins were simulated with charged termini. The tripeptides were placed in a periodic rectangular box of dimensions 5 × 5 × 5 nm³ with approximately 4000 CHARMM TIP3P (59,60) water molecules in the AA simulations and 950 polarizable water particles in the CG simulations. (61) The water-soluble protein cardiotoxin V was placed in a periodic rectangular box of 7.9 × 7.9 × 7.9 nm³ and filled with 16500 CHARMM TIP3P water molecules in the AA simulations. In the CG simulations, the protein was placed inside a periodic rectangular box of 5.7 × 5.7 × 5.7 nm³ and filled with 1800 polarizable water particles. The larger water-soluble protein HEWL was placed in a periodic rectangular box of 8.9 × 8.9 × 8.9 nm³ and filled with 23000 CHARMM TIP3P water molecules in the AA simulations and 5400 polarizable water particles in the CG simulations. Na⁺ and Cl^– ions were added to all systems at 0.15 M concentration to neutralize the protein systems. The turkey ovomucoid inhibitor protein was placed in a box of 4.9 × 4.9 × 4.9 nm³ with 3086 SPC (62) water molecules. The GLIC protein was embedded into a bilayer membrane containing 498 phosphatidylcholine (POPC) lipids, placed in a box of 14.0 × 14.0 × 15.9 nm³, and filled with 66494 CHARMM TIP3P waters, 58 Na⁺, and 123 Cl^– ions. The system contained 292135 atoms in total. The simulation of this system was performed with the GROMACS 2021.4 release as reference. The GLIC benchmarks were run with default settings on an Intel i9–7920X 12-core CPU and an Nvidia RTX 2080 Ti GPU. All input configurations are provided as the Supporting Information.

In the AA simulations, Coulomb interactions were modeled with the smooth PME method with a real-space cutoff of 1.2 nm and a grid spacing of 0.14 nm, (44,45) while Lennard-Jones interactions were smoothly switched to zero in a range from 1.0 to 1.2 nm. In the CG simulations, Coulomb interactions were modeled by a Reaction Field potential with a 1.1 nm cutoff, ϵ_r = 2.5, and ϵ_RF = ∞, while Lennard-Jones interactions were truncated at 1.1 nm. (63) To keep the temperature constant at 300 K, we used the v-rescale thermostat (64) with time constants of 0.5 and 1.0 ps^–1 for AA and CG simulations, respectively. The pressure was kept constant at 1 bar with the Parrinello–Rahman barostat (65) with relaxation times of 2.0 and 12.0 ps for AA and CG simulations, respectively. A leapfrog integrator was used with an integration step of 2 and 20 fs for AA and CG simulations, respectively. In the AA simulations, the LINCS algorithm (66) was used to constrain h-bond lengths of the solutes, while the SETTLE (67) algorithm was used to constrain internal degrees of freedom of the water molecules. Prior to the constant pH MD simulations, the potential energy of each system was minimized using the steepest descent method, followed by 1 ns of equilibration.

Constant pH MD Simulation Setups

In the atomistic simulations, the multisite representation was used to model the protonation states of titratable residues. Two λ-coordinates were introduced to model the two forms of the carboxylic acid side chain in Asp and Glu, while three coordinates were used to describe the three protonation states of the imidazole side chain in His. In the CG simulations, the single-site representation was used, in which the A and B states represent the protonated and deprotonated states of the titratable beads. Because, in contrast to AA force fields, there is no distinction between the two neutral forms of the His side chain in the Martini force field, the single-site description for HIS suffices in the CG simulations. In both atomistic and coarse-grained simulations, the transformations between the different protonation states were achieved by changing the charges of the ionizable groups. The Lennard-Jones and bonded terms (bonds, angles, and torsions) were kept in the protonated and deprotonated states in AA and CG simulations, respectively. We show in Figure S3 that the contribution of these terms is sufficiently small to be neglected without significant error. We note, however, that these terms can be made λ-dependent as well, but this is beyond the scope of the current work since the efforts to implement this are high.

The mass of the λ-particles was set to 5 atomic units, and their temperature was maintained at 300 K by using a separate v-rescale thermostat for the λ-coordinates with a time constant of 2.0 ps^–1. For all λ-coordinates the biasing potential V_i^bias(λ_i) was defined by equation S1 in the Supporting Information. The barrier height of the double-well potential was set to 5.0 and 7.5 kJ/mol for AA and CG simulations, respectively. The parameters for the double-well potential and the pH-dependent potential (eq 4) are provided in Table S1.

For the tripeptides, we calculated five independent CpHMD trajectories of 20 ns each at 13 pH values, ranging from 1.0 to 7.0 for the peptides with Glu and Asp, and from 4.0 to 10.0 for the peptides with His. For the cardiotoxin V protein (three Asp and one His titratable residues), we performed five independent CpHMD simulations of 50 ns at 15 pH values between 1.0 to 8.0. For the HEWL protein (seven Asp, two Glu, and one His titratable residues), we performed five independent CpHMD simulations of 75 ns at 21 pH values between −1.0 to 9.0. The values of the λ-coordinates were written to the output file with a frequency of 1 ps^–1.

Reference States and Force Field Correction Potentials

The constant pH simulations of the aforementioned systems require reference states for Asp, Glu, and His, in which the proton affinity (pK_a) is known from the experiment. The measured and calculated (force field) deprotonation free energies of these reference states were used to include the effect of the pH bath, V^pH(λ), as well as the effects of the breaking and forming of chemical bonds in the simulation, i.e., V^MM in eq 2. The measured reference pK_a values used in this work are included in Table 1. Note that the experimental values were obtained for pentapeptides, while tripeptides were used for computing V^MM. This however did not affect the results, as shown in Figure S2.

Table 1. pK_a Values Obtained from Titration Simulations^a

	tripeptide simulations (69,70)
	pK_a values
amino acid	CHARMM36	MARTINI	exp.
Asp	3.61 ± 0.03	3.69 ± 0.02	3.65
Glu	4.26 ± 0.04	4.30 ± 0.03	4.25
His macroscopic	6.34 ± 0.08	6.40 ± 0.03	6.42
His HSD	6.56 ± 0.06		6.53
His HSE	6.90 ± 0.05		6.94

	simulation of cardiotoxin V (71,72)
	pK_a values
amino acid	CHARMM36	MARTINI	exp.
His-4	5.14 ± 0.16	4.54 ± 0.09	5.5
Glu-17	4.08 ± 0.08	4.36 ± 0.04	4
Asp-42	4.02 ± 0.10	4.30 ± 0.05	3.2
Asp-59	2.41 ± 0.07	1.45 ± 0.03	<2
	r = 0.96	r = 0.80
	MSE = 0.24	MSE = 0.64
	RMSE = 0.49	RMSE = 0.80

	simulation of HEWL (73)
	pK_a values
amino acid	CHARMM36	MARTINI	exp.
Glu-7	2.82 ± 0.07	4.86 ± 0.05	2.6
His-15	4.84 ± 0.05	5.42 ± 0.05	5.5
Asp-18	3.35 ± 0.05	3.31 ± 0.03	2.8
Glu-35	7.64 ± 0.13	6.36 ± 0.05	6.1
Asp-48	0.99 ± 0.07	3.36 ± 0.05	1.4
Asp-52	5.69 ± 0.12	7.18 ± 0.10	3.6
Asp-66	1.70 ± 0.10	5.22 ± 0.05	1.2
Asp-87	1.73 ± 0.03	3.47 ± 0.05	2.2
Asp-101	5.43 ± 0.11	4.20 ± 0.06	4.5
Asp-119	2.77 ± 0.05	3.80 ± 0.05	3.5
	r = 0.90	r = 0.49
	MSE = 0.96	MSE = 4.01
	RMSE = 0.98	RMSE = 2.00

^{^a}

The reference pK_a values for tripeptides are given in the last column that contain the experimental pK_a values. The values for Asp and Glu are taken from ref (69), while the microscopic and macroscopic pK_a values for His are taken from ref (70). Experimental pK_a values for cardiotoxin V are from refs (71and72) and for HEWL from ref (73). For both proteins Pearson correlation (r), MSE and RMSE errors are provided.

Thermodynamic integration was used to compute the reference free energies as follows: the partial charges in tripeptide systems representing the reference states of Glu, Asp, and His were linearly interpolated between λ = −0.1 and λ = 1.1 with a step of 0.05 under the constraint λ₁ + λ₂ = 1 for Glu and Asp, while for His, the constraint was λ₁ + λ₂ + λ₃ = 1. For each set of λ values, called a grid point, a 10 ns MD simulation was performed. The ∂V/∂λ_i values were saved every ps, which is approximately equal to the autocorrelation times for the λ-coordinates. The total charge of the system was kept neutral by simultaneously changing the charge of a single buffer particle, as discussed below. The ∂V/∂λ_i values were averaged over the last 9 ns of the trajectories. To obtain an analytical expression for V^MM, a fifth-order polynomial was fitted to these averages for Asp and Glu, while an eighth-order polynomial was fitted for His, taking into account possible linear dependencies of the coefficients (see the Supporting Information). Fitting errors were below 0.5 kJ/mol for Asp and Glu and below 1 kJ/mol for His, which are of similar magnitude as the statistical accuracy of the derivatives.

Buffer Particles

Dynamically changing partial charges can affect the total charge of the simulation unit cell, which can lead to artifacts, as documented for instance in Hub et al. for Ewald-based methods. (68) To avoid such artifacts, it is essential to keep the total charge of the unit cell constant. Two approaches have been proposed: (i) direct coupling between each titratable residue and a water, (27) or ion, (25) and (ii) titratable buffers that collectively compensate for changes in charge of all titratable residues. (42)

Here, we follow the latter approach, but with several improvements for all-atom simulations. First, to avoid restraints, which were needed to minimize interactions between the buffers and the titrable sites in previous work, (42) we introduced buffer particles with both small LJ radius and small partial charges of maximal |0.5|e, such that they do not disturb the hydrogen bond network, nor interact too strongly with the titratable sites or other buffers. Second, to also prevent strong interactions with hydrophobic regions in the system, the C⁽⁶⁾ dispersion parameter with anything other than water was set to zero, including the other buffers. The latter also avoids the clustering of buffers during the simulation. Thus, the buffer particles have an σ of 0.25 nm and an ϵ of 4 kJ/mol. Further details on buffer parametrization are provided in the accompanying paper. (55) In coarse-grained simulations, standard Na⁺ ions were used as buffer particles.

As in Donnini et al., (42) the buffers were collectively coupled to the titratable sites in the system via a charge constraint. The charges of all buffers were thus simultaneously interpolated between −0.5e and 0.5e, keeping the simulation box neutral. For all peptide simulations, 10 such buffers were introduced into the system, while 20 and 50 buffers were added to the simulation boxes with cardiotoxin V and HEWL proteins (systems 4–5), respectively, in both AA and CG models. 185 buffer particles were used in GLIC simulations.

Analysis of the Constant pH Trajectories

To estimate the pK_a values of titratable groups from multiple simulations at various pH values, we computed the average fraction of deprotonated frames (S^deprot) over all replicas. For a group with a single titratable site, this average was obtained as

S^{deprot} (pH) = \frac{N^{deprot}}{N^{prot} + N^{deprot}}

(18)

where N^prot and N^deprot are the total number of frames in which the site is protonated and deprotonated, respectively. For titratable sites modeled in the single-site representation, we considered it protonated if λ is below 0.2 and deprotonated if λ is above 0.8. For sites that are described with the multisite description, we considered a state protonated if the λ associated with the protonated form of the residue is above 0.8 and deprotonated if the λ associated with the deprotonated form of the residue is above 0.8.

To estimate the macroscopic pK_a values of histidine, which contains two titratable sites N_ϵ and N_δ, we calculated for each pH value the average fraction of frames in which the residue is deprotonated at either of the two sites:

S_{macro}^{deprot} (pH) = 1 - \frac{N^{λ_{p}}}{N^{λ_{p}} + N^{λ_{ϵ}} + N^{λ_{δ}}}

(19)

where N^λ_p, N^λ_ϵ, and N^λ_δ are the numbers of frames in which λ_p > 0.8, λ_ϵ > 0.8, and λ_δ > 0.8 (Figure 2). To estimate the microscopic pK_a values for the two sites of His, we calculated for each site the average fraction of frames in which that site was deprotonated:

S_{micro (δ / ϵ)}^{deprot} (pH) = \frac{N^{λ_{δ} / λ_{ϵ}}}{N^{λ_{p}} + N^{λ_{ϵ} / λ_{δ}}}

(20)

Errors were estimated from the standard error of the mean for the different replicas.

The averaged fractions at each pH value were fitted to the Henderson–Hasselbalch equation:

S^{deprot} = \frac{1}{10^{(p K_{a} - pH)} + 1}

(21)

which yielded the pK_a values as fitting parameters. The error in the pK_a was estimated from the 95% confidence interval for the nonlinear least-squares fit to the average (S^deprot) values.

Results and Discussion

Click to copy section linkSection link copied!

Here we discuss the results obtained with our new implementation of constant pH into the fork of the GROMACS 2021 release. (49) While here our focus is on the validity and performance of the constant pH MD implementation, the convergence of the conformational and λ degrees of freedom are investigated systematically in the accompanying paper. (55)

Titration of Single Amino Acids

In Figure 3, we show the titration curves for AlaAspAla, AlaGluAla, and AlaHisAla tripeptides, obtained from simulations with the modified all-atom CHARMM36 (55) and coarse-grained Martini 2.0 force fields. (41) Fitting the deprotonated fractions as a function of pH value to the Henderson–Hasselbalch equation (dashed lines in Figure 3) yields pK_a values for the tripeptides that are within 0.1 pK_a units from the reference values. Comparing the titration curves obtained with the Martini 2.0 force field in our implementation to those computed with the constant pH approach developed explicitly for this coarse-grained model, (34) our results suggest a much better agreement with the experiment than the latter. We attribute this difference to the more sophisticated explicit treatment of proton-like particles in the Martini constant pH approach. The rather good agreement between the titration curves obtained for both force fields on the one hand and the experiment on the other hand suggests that our implementation has little to no dependency on the force field, in line with the GROMACS philosophy of supporting a wide range of popular force fields.

Figure 3. Titrations of tripeptide amino acids (Glu, Asp, and His) in water. The top and bottom rows show titrations performed with the modified AA CHARMM36 (55) and CG Martini (41) force fields, respectively. In all simulations, neutrality was maintained by including ten buffer particles in combination with the charge constraint. Dots show the fraction of frames in which the residue is deprotonated, and the dashed lines represent the fits to the Henderson–Hasselbalch equation. For His, the blue color represents the macroscopic pK_a, while yellow and red represent the microscopic pK_a values for HSD (proton on N_δ) and HSE (proton on N_ϵ), respectively. In the Martini 2.0 model, HSD and HSE are indistinguishable and hence only the macroscopic titration curve is shown. Errors of S^deprot were estimated from the standard error of the mean for the different replicas. From the fits, the pK_a values were estimated and listed in Table 1.

Titration of Proteins

The titration curves of cardiotoxin V and HEWL proteins are shown in Figures 4 and 5, respectively. The pK_a values obtained from fitting the Henderson–Hasselbalch equation to the degree of deprotonation in the all-atom simulations of both proteins with the CHARMM36 force field, listed in Table 1, are in good agreement with previous constant pH MD simulations (30,74) and in reasonable agreement with experimental estimates from NMR spectroscopy [Pearson correlation coefficient (r) 0.96 and 0.9, RMSE 0.49 and 0.98 for cardiotoxin V and HEWL, respectively]. (71−73) The pK_a values and S^deprot are converged in 50 ns, as discussed in the Supporting Information section 5 (Figures S6–S9). Analysis of the RMS deviation of the backbone and of the RMS fluctuation of the residues, plotted in Figures S10–S13, suggest no major influence of the pH on the structural stability of these proteins. The titration correlates with solvent exposure, which contributes to the stabilization of the charged protonation state (see Supporting Information section 7) (Figures S14–S16). The trends in the pK_a shifts are well reproduced, including the downshift of Asp-59 in cardiotoxin V, and, with the exception of Glu-35 and Asp-52 in HEWL, the deviations are below 1 pK_a unit. We note that also in previous constant pH simulations with the CHARMM force field, (22,30) similar deviations were found for these two residues (see Figure S5). This suggests that the origin of the discrepancy might lie beyond the implementation, and could be due to either a lack of sampling or systematic shortcomings in the force field, as was discussed in Huang et al. (30) For detailed insights into the structural origins of these pK_a shifts, we refer the reader to the paper of Swails and Roitberg. (19)

Figure 4. Titration curves of the cardiotoxin V protein obtained from constant pH MD simulations with the CHARMM36 (top) and Martini 2.0 force fields (bottom). For each of the four titratable residues in this protein, the dots show the fraction of frames in which the residue is deprotonated. Errors of S^deprot were estimated from the standard error of the mean for the different replicas. The lines show the best fits to the Henderson–Hasselbalch equation. The pK_a values for each titratable residue were estimated from these fits and listed in Table 1. The right panel shows the protein structure with the four titratable residues highlighted in stick representation.

Figure 5. Titration curves of the HEWL protein obtained from constant pH MD simulations with the CHARMM36 (top) and Martini 2.0 force fields (bottom). For each of the ten titratable residues, the dots show the fraction of frames in which that residue is deprotonated. Errors of S^deprot were estimated from the standard error of the mean for the different replicas. The lines show the best fit to the Henderson-Hasselbach equation. The pK_a values for each titratable residue were estimated from these fits and listed in Table 1. The right panel shows the protein structure with the ten titratable residues highlighted in stick representation.

The pK_a values estimated from the Martini 2.0 force field simulations of these proteins do not agree as well with the experiment as those derived from the all-atom simulations (Pearson correlation coefficient (r) 0.8 and 0.49; RMSE 0.8 and 2.0 for cardiotoxin V and HEWL, respectively). We speculate that the larger deviation of the pK_a’s in the coarse-grained constant pH simulations could be due to the lower accuracy of the electrostatic interactions. Although we still consider the results obtained with the Martini simulations reasonable, in particular for the peptides, the discrepancies for the titratable residues in proteins suggest that additional parametrization efforts may be required to systematically improve the force field for constant pH MD simulations based on λ-dynamics. Such improvements would be particularly worthwhile considering coarse grained simulations pave the way to perform MD simulations of complete organelles, (75) in which many processes have a strong pH dependence.

Efficiency of the Implementation

To demonstrate that linear interpolation of charges is computationally more efficient than the linear interpolation of the potential energy functions for systems with many titratable sites, we investigated how the computational cost of a simulation scales with the number of titratable sites in the system for both approaches. Because we have implemented the interpolation of the charges rather than potential energy functions into the fork of GROMACS 2021 release, whereas the potential energy function interpolation was implemented in a fork of GROMACS 3.3 release, we compare the relative performances of both codes for an increasing number of titratable sites in the system. We define the relative performance as the ratio between the average number of integration steps per time unit for a simulation with constant pH on the one hand and the average number of integration steps per time unit for a normal simulation without constant pH on the other hand.

Figure 6A shows that the relative performance of constant pH simulations with charge interpolation does not decrease when the number of titratable sites included in the simulation increases. Most of the 30–40% drop in performance compared to a normal MD simulation with the same version of GROMACS is caused by the additional calculations and reductions in the nonbonded pair-interaction kernels that are required to obtain the real-space part of the electrostatic potential (Φ(R_i, λ) in eq 17).

Figure 6. (A) Relative performance of interpolating potentials in a previous implementation of CpHMD into a fork of GROMACS 3.3 release (blue) and of charge interpolation in our new implementation (red) as a function of the number of titratable sites. The simulations were performed for the turkey ovomucoid inhibitor protein (PDB ID: 2GKR (53)), shown in the inset, where the titratable sites are highlighted in stick representation. (B) Comparison of the performance between CPU-only and CPU+GPU implementations for the ligand-gated ion channel GLIC (PDB ID: 4HFI (52)) with 185 titratable sites. In total, the GLIC system contained 292135 atoms.

In contrast, the relative performance of constant pH simulations based on the linear interpolation of potential energy functions decreases with the number of titratable sites in the system. This comparison thus demonstrates that by replacing linear interpolation of potentials with linear interpolation of partial charges, we have overcome the major bottleneck in the earlier constant pH implementation in the fork of GROMACS 3.3 release and paved the way toward simulations of large biomolecular systems at constant pH.

An example of such a large system is the proton-gated ion channel GLIC, a membrane protein with 185 titratable residues. Figure 6B shows the performance of the new implementation for this large system when running the simulation on CPU and on a combination of CPU and GPU. While the computational overhead is somewhat larger when using a GPU in addition to a CPU, the overall performance still improves significantly when adding a GPU. We have also implemented a parallel version using MPI. For the GLIC system, the code scales up to 256 cores, on the Mahti supercomputer at CSC, with a performance of 42 ns/day, compared with 61 ns/day without constant pH.

Conclusions

Click to copy section linkSection link copied!

We have presented and validated a new implementation of λ-dynamics based constant pH molecular dynamics in the GROMACS software. Our implementation combines several developments in this field into a single MD program, including the multisite representation of titratable groups, (46) charge interpolation, (21) Particle Mesh Ewald electrostatics, (30) and charge constraints. (42) Test calculations on amino acids and proteins suggest that the new implementation is efficient, accurate, and agnostic to force fields. Combined with user-friendly parametrization protocols, presented in the accompanying paper, (55) we expect that this implementation will pave the way toward routinely including the effect of pH in biomolecular MD simulations. (76)

Supporting Information

Click to copy section linkSection link copied!

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.2c00516.

The input files and parameters used for MD simulations in the presented work (ZIP)
A Mathematica notebook with instructions and routines to fit V^MM (ZIP)
Description of λ-potentials, the effect of neglecting the interpolation of Lennard-Jones interactions, titration results for Asp and Glu within the single-site representation, a comparison of pK_a values for HEWL obtained with various λ-dynamics-based constant pH methods, demonstration that charge interpolation requires a single evaluation of the electrostatic potential for both single- and multisite representations (PDF)

Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Author Information

Click to copy section linkSection link copied!

Corresponding Authors
- Noora Aho - Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014Jyväskylä, Finland; Email: [email protected]
- Pavel Buslaev - Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014Jyväskylä, Finland; https://orcid.org/0000-0003-2031-4691; Email: [email protected]
- Gerrit Groenhof - Nanoscience Center and Department of Chemistry, University of Jyväskylä, 40014Jyväskylä, Finland; https://orcid.org/0000-0001-8148-5334; Email: [email protected]
- Berk Hess - Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44Stockholm, Sweden; Email: [email protected]
Authors
- Anton Jansen - Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44Stockholm, Sweden
- Paul Bauer - Department of Applied Physics and Swedish e-Science Research Center, Science for Life Laboratory, KTH Royal Institute of Technology, 100 44Stockholm, Sweden; https://orcid.org/0000-0002-2268-0065
Author Contributions
N.A. and P.B. contributed equally.
Notes
The authors declare no competing financial interest.
The fork of GROMACS 2021 with constant pH implemented as described here, is available for download free of charge from https://gitlab.com/gromacs-constantph/constantph. In addition to the source code, also instructions on how to set up and perform MD simulations are available

Acknowledgments

Click to copy section linkSection link copied!

This research was supported by the Swedish Research Council (Grant 2019-04477), Academy of Finland (Grants 311031 and 332743), and the BioExcel CoE (Grant H2020-INFRAEDI-02-2018-823830). The simulations were performed on resources provided by the CSC ─ IT Center for Science, Finland, and the Swedish National Infrastructure for Computing (SNIC 2021/1-38). We also thank Erik Lindahl, Helmut Grubmüller, Dmitry Morozov, Serena Donnini, and Plamen Dobrev for their support during the project.

Appendix: Constraint Algorithm

Click to copy section linkSection link copied!

We use constraints to restrict sampling to the correct protonation states in the multisite representation as well as to maintain a neutral charge of the simulation box. Both multisite and charge constraints keep the linear combination of a subset of λ-coordinates constant and are applied simultaneously. Thus, there are N_c constraint equations

σ^{k} (λ) = \sum_{i = 1}^{N_{sites}} α_{i}^{k} λ_{i} = C^{k}

(22)

for k ≤ N_c. Here, λ is the vector of all λ_i-coordinates, and C^k is the value of constraint k, which can be zero. If σ^k(λ) is a multisite constraint, α_i^k = 1 for λ_i-coordinates that represent one of the protonation states of a residue, while α_i^k = 0 for all other λ_i-coordinates. If σ^k(λ) is a constraint for keeping the overall charge constant, α_i^k = ∑_j^N_atomsⁱq_j,i^B – q_j,i^A, with N_atomsⁱ the number of atoms whose charges change as a function of λ_i.

During a leapfrog integration step, all λ_i-coordinates are first propagated without constraints to their unconstrained new values λ_i^u, which do not fulfill the constraints in eq 22. To obtain the constrained λ_i^c values, we first connect the constrained and unconstrained λ_i values using the definition of σ^k(λ):

σ^{k} (λ^{c}) = σ^{k} (λ^{u}) + \sum_{i = 1}^{N_{sites}} α_{i}^{k} [λ_{i}^{c} - λ_{i}^{u}]

(23)

Because the unconstrained and constrained λ_i-coordinates are also connected by the constraint forces

G_{i}^{k} = - ζ^{k} \frac{\partial σ^{k}}{\partial λ_{i}}

, we have in addition that

λ_{i}^{c} = λ_{i}^{u} + \sum_{k = 1}^{N_{c}} G_{i}^{k} \frac{Δ t^{2}}{m_{i}} = λ_{i}^{u} - \sum_{k = 1}^{N_{c}} ζ^{k} \frac{\partial σ^{k}}{\partial λ_{i}} \frac{Δ t^{2}}{m_{i}} = λ_{i}^{u} - \sum_{k = 1}^{N_{c}} ζ^{k} α_{i}^{k} \frac{Δ t^{2}}{m_{i}}

(24)

where ζ^k is the Lagrange multiplier for constraint k, m_i the fictitious mass of λ_i, and Δt the integration time step. Substituting eq 24 in 23 yields

σ^{k} (λ^{c}) = σ^{k} (λ^{u}) - \sum_{i = 1}^{N_{sites}} α_{i}^{k} \sum_{l = 1}^{N_{c}} ζ^{l} α_{i}^{l} \frac{Δ t^{2}}{m_{i}}

(25)

which after rearranging can be expressed as

Δ σ^{k} = σ^{k} (λ^{u}) - σ^{k} (λ^{c}) = \sum_{l = 1}^{N_{c}} ζ^{l} \sum_{i = 1}^{N_{sites}} α_{i}^{k} α_{i}^{l} \frac{Δ t^{2}}{m_{i}}

(26)

The last expression can be rewritten in matrix form

Δ σ = {(Δ σ^{1}, ..., Δ σ^{N_{c}})}^{T} = A ζ

(27)

where ζ = (ζ¹,···,ζ^N_c)^T and

A = (\sum_{i = 1}^{N_{sites}} α_{i}^{k} α_{i}^{l} \frac{Δ t^{2}}{m_{i}})

(28)

Because the α_i^k coefficients remain the same, matrix A is computed once at the start of the simulation. At each step, the elements Δσ^k are evaluated as the difference between the σ^k(λ^u) and σ^k(λ^c) = C^k:

Δ σ^{k} = σ^{k} (λ^{u}) - C^{k}

(29)

The Lagrange multipliers ζ^k are then obtained from eq 27 and used to correct the unconstrained λ_i values (eq 24).

References

Click to copy section linkSection link copied!

This article references 76 other publications.

1
Hollingsworth, S. A.; Dror, R. O. Molecular Dynamics Simulation for All. Neuron 2018, 99, 1129– 1143, DOI: 10.1016/j.neuron.2018.08.011
Google Scholar
1
Molecular Dynamics Simulation for All
Hollingsworth, Scott A.; Dror, Ron O.
Neuron (2018), 99 (6), 1129-1143CODEN: NERNET; ISSN:0896-6273. (Cell Press)
A review. The impact of mol. dynamics (MD) simulations in mol. biol. and drug discovery has expanded dramatically in recent years. These simulations capture the behavior of proteins and other biomols. in full at. detail and at very fine temporal resoln. Major improvements in simulation speed, accuracy, and accessibility, together with the proliferation of exptl. structural data, have increased the appeal of biomol. simulation to experimentalists-a trend particularly noticeable in, although certainly not limited to, neuroscience. Simulations have proven valuable in deciphering functional mechanisms of proteins and other biomols., in uncovering the structural basis for disease, and in the design and optimization of small mols., peptides, and proteins. Here we describe, in practical terms, the types of information MD simulations can provide and the ways in which they typically motivate further exptl. work.
2
Groenhof, G.; Modi, V.; Morozov, D. Observe while it happens: catching photoactive proteins in the act with non-adiabatic molecular dynamics simulations. Curr. Opin. Struct. Biol. 2020, 61, 106– 112, DOI: 10.1016/j.sbi.2019.12.013
Google Scholar
2
Observe while it happens: catching photoactive proteins in the act with non-adiabatic molecular dynamics simulations
Groenhof, Gerrit; Modi, Vaibhav; Morozov, Dmitry
Current Opinion in Structural Biology (2020), 61 (), 106-112CODEN: COSBEF; ISSN:0959-440X. (Elsevier Ltd.)
A review. Organisms use photo-receptors to react to light. The first step is usually the absorption of a photon by a prosthetic group embedded inside the photo-receptor, often a conjugated chromophore. The electronic changes in the chromophore induced by photo-absorption can trigger a cascade of structural or chem. transformations that culminate into a response to light. Understanding how these proteins have evolved to mediate their activation process has remained challenging because the required time and spacial resolns. are notoriously difficult to achieve exptl. Therefore, mechanistic insights into photoreceptor activation have been predominantly obtained with computer simulations. Here we briefly outline the challenges assocd. with such computations and review the progress made in this field.
3
Warshel, A.; Sussman, F.; King, G. Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process. Biochemistry 1986, 25, 8368– 8372, DOI: 10.1021/bi00374a006
Google Scholar
3
Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process
Warshel, A.; Sussman, F.; King, G.
Biochemistry (1986), 25 (26), 8368-72CODEN: BICHAW; ISSN:0006-2960.
Evaluation of the free energy of ionization of acidic groups in proteins may be used as a powerful and general test case for detg. the reliability of calcns. of electrostatic energies in macromols. This test case was addressed by using an adiabatic charging process that evaluates the changes in free energies assocd. with ionizing the acidic groups of aspartate-3 and glutamate-7 in bovine pancreatic trypsin inhibitor and aspartic acid in soln. The results of these free energy calcns. are very encouraging; the error range is ∼1 kcal/mol for these values, which are ∼70 kcal/mol. Thus, the stage of obtaining quant. results in modeling the energetics of solvated proteins is finally being approached.
4
Alexov, E.; Mehler, E. L.; Baker, N.; Baptista, A. M.; 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 pKa values in proteins. Proteins: Struct., Funct., Bioinf. 2011, 79, 3260– 3275, DOI: 10.1002/prot.23189
Google Scholar
4
Progress in the prediction of pKa values in proteins
Alexov, Emil; Mehler, Ernest L.; Baker, Nathan; Baptista, Antonio M.; Huang, Yong; Milletti, Francesca; Erik Nielsen, Jens; Farrell, Damien; Carstensen, Tommy; Olsson, Mats H. M.; Shen, Jana K.; Warwicker, Jim; Williams, Sarah; Word, J. Michael
Proteins: Structure, Function, and Bioinformatics (2011), 79 (12), 3260-3275CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
A review. The pKa-cooperative aims to provide a forum for exptl. and theor. researchers interested in protein pKa values and protein electrostatics in general. The first round of the pKa-cooperative, which challenged computational labs to carry out blind predictions against pKas exptl. detd. in the lab. of Bertrand Garcia-Moreno, was completed and results discussed at the Telluride meeting (July 6-10, 2009). This article serves as an introduction to the reports submitted by the blind prediction participants that will be published in a special issue of PROTEINS: Structure, Function and Bioinformatics. Here, the authors briefly outline existing approaches for pKa calcns., emphasizing methods that were used by the participants in calcg. the blind pKa values in the first round of the cooperative. The authors then point out some of the difficulties encountered by the participating groups in making their blind predictions, and finally try to provide some insights for future developments aimed at improving the accuracy of pKa calcns. Proteins 2011; © 2011 Wiley-Liss, Inc.
5
Chen, W.; Morrow, B. H.; Shi, C.; Shen, J. K. Recent development and application of constant pH molecular dynamics. Mol. Simul. 2014, 40, 830– 838, DOI: 10.1080/08927022.2014.907492
Google Scholar
5
Recent development and application of constant pH molecular dynamics
Chen, Wei; Morrow, Brian H.; Shi, Chuanyin; Shen, Jana K.
Molecular Simulation (2014), 40 (10-11), 830-838CODEN: MOSIEA; ISSN:0892-7022. (Taylor & Francis Ltd.)
A review. Soln. pH is a crit. environmental factor for chem. and biol. processes. Over the last decade, significant efforts have been made in the development of const. pH mol. dynamics (pHMD) techniques for gaining detailed insights into pH-coupled dynamical phenomena. In this article, we review the advancement of this field in the past five years, placing a special emphasis on the development of the all-atom continuous pHMD technique. We discuss various applications, including the prediction of pKa shifts for proteins, nucleic acids and surfactant assemblies, elucidation of pH-dependent population shifts, protein-protein and protein-RNA binding, as well as the mechanisms of pH-dependent self-assembly and phase transitions of surfactants and peptides. We also discuss future directions for the further improvement of the pHMD techniques.
6
Zeng, X.; Mukhopadhyay, S.; Brooks, C. L., III Residue-level resolution of alphavirus envelope protein interactions in pH-dependent fusion. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 2034– 2039, DOI: 10.1073/pnas.1414190112
Google Scholar
6
Residue-level resolution of alphavirus envelope protein interactions in pH-dependent fusion
Zeng, Xiancheng; Mukhopadhyay, Suchetana; Brooks, Charles L., III
Proceedings of the National Academy of Sciences of the United States of America (2015), 112 (7), 2034-2039CODEN: PNASA6; ISSN:0027-8424. (National Academy of Sciences)
Alphavirus envelope proteins, organized as trimers of E2-E1 heterodimers on the surface of the pathogenic alphavirus, mediate the low pH-triggered fusion of viral and endosomal membranes in human cells. The lack of specific treatment for alphaviral infections motivates our exploration of potential antiviral approaches by inhibiting one or more fusion steps in the common endocytic viral entry pathway. In this work, we performed const. pH mol. dynamics based on an at. model of the alphavirus envelope with icosahedral symmetry. We have identified pH-sensitive residues that cause the largest shifts in thermodn. driving forces under neutral and acidic pH conditions for various fusion steps. A series of conserved interdomain His residues is identified to be responsible for the pH-dependent conformational changes in the fusion process, and ligand binding sites in their vicinity are anticipated to be potential drug targets aimed at inhibiting viral infections.
7
Law, S. M.; Zhang, B. W.; Brooks, C. L., III pH-sensitive residues in the p19 RNA silencing suppressor protein from carnation Italian ringspot virus affect siRNA binding stability. Protein Sci. 2013, 22, 595– 604, DOI: 10.1002/pro.2243
Google Scholar
7
pH-sensitive residues in the p19 RNA silencing suppressor protein from carnation Italian ringspot virus affect siRNA binding stability
Law, Sean M.; Zhang, Bin W.; Brooks, Charles L., III
Protein Science (2013), 22 (5), 595-604CODEN: PRCIEI; ISSN:1469-896X. (Wiley-Blackwell)
Tombusviruses, such as Carnation Italian ringspot virus (CIRV), encode a protein homodimer called p19 that is capable of suppressing RNA silencing in their infected hosts by binding to and sequestering short-interfering RNA (siRNA) away from the RNA silencing pathway. P19 binding stability has been shown to be sensitive to changes in pH but the specific amino acid residues involved have remained unclear. Using const. pH mol. dynamics simulations, we have identified key pH-dependent residues that affect CIRV p19-siRNA binding stability at various pH ranges based on calcd. changes in the free energy contribution from each titratable residue. At high pH, the deprotonation of Lys60, Lys67, Lys71, and Cys134 has the largest effect on the binding stability. Similarly, deprotonation of several acidic residues (Asp9, Glu12, Asp20, Glu35, and/or Glu41) at low pH results in a decrease in binding stability. At neutral pH, residues Glu17 and His132 provide a small increase in the binding stability and we find that the optimal pH range for siRNA binding is between 7.0 and 10.0. Overall, our findings further inform recent expts. and are in excellent agreement with data on the pH-dependent binding profile.
8
Ellis, C. R.; Shen, J. pH-dependent population shift regulates BACE1 activity and inhibition. J. Am. Chem. Soc. 2015, 137, 9543– 9546, DOI: 10.1021/jacs.5b05891
Google Scholar
8
pH-Dependent Population Shift Regulates BACE1 Activity and Inhibition
Ellis, Christopher R.; Shen, Jana
Journal of the American Chemical Society (2015), 137 (30), 9543-9546CODEN: JACSAT; ISSN:0002-7863. (American Chemical Society)
BACE1, a major therapeutic target for treatment of Alzheimer's disease, functions within a narrow pH range. Despite tremendous effort and progress in the development of BACE1 inhibitors, details of the underlying pH-dependent regulatory mechanism remain unclear. Here the authors elucidate the pH-dependent conformational mechanism that regulates BACE1 activity using continuous const.-pH mol. dynamics (MD). The simulations reveal that BACE1 mainly occupies three conformational states and that the relative populations of the states shift according to pH. At intermediate pH, when the catalytic dyad is monoprotonated, a binding-competent state is highly populated, while at low and high pH a Tyr-inhibited state is dominant. The authors' data provide strong evidence supporting conformational selection as a major mechanism for substrate and peptide-inhibitor binding. These new insights, while consistent with expt., greatly extend the knowledge of BACE1 and have implications for further optimization of inhibitors and understanding potential side effects of targeting BACE1. Finally, the work highlights the importance of properly modeling protonation states in MD simulations.
9
Kim, M. O.; Blachly, P. G.; McCammon, J. A. Conformational dynamics and binding free energies of inhibitors of BACE-1: from the perspective of protonation equilibria. PLoS computational biology 2015, 11, e1004341, DOI: 10.1371/journal.pcbi.1004341
Google Scholar
9
Conformational dynamics and binding free energies of inhibitors of BACE-1: from the perspective of protonation equilibria
Kim, M. Olivia; Blachly, Patrick G.; McCammon, J. Andrew
PLoS Computational Biology (2015), 11 (10), e1004341/1-e1004341/28CODEN: PCBLBG; ISSN:1553-7358. (Public Library of Science)
BACE-1 is the β-secretase responsible for the initial amyloidogenesis in Alzheimer's disease, catalyzing hydrolytic cleavage of substrate in a pH-sensitive manner. The catalytic mechanism of BACE-1 requires water-mediated proton transfer from aspartyl dyad to the substrate, as well as structural flexibility in the flap region. Thus, the coupling of protonation and conformational equil. is essential to a full in silico characterization of BACE-1. In this work, we perform const. pH replica exchange mol. dynamics simulations on both apo BACE-1 and five BACE-1-inhibitor complexes to examine the effect of pH on dynamics and inhibitor binding properties of BACE-1. In our simulations, we find that soln. pH controls the conformational flexibility of apo BACE-1, whereas bound inhibitors largely limit the motions of the holo enzyme at all levels of pH. The microscopic pKa values of titratable residues in BACE-1 including its aspartyl dyad are computed and compared between apo and inhibitor-bound states. Changes in protonation between the apo and holo forms suggest a thermodn. linkage between binding of inhibitors and protons localized at the dyad. Utilizing our recently developed computational protocol applying the binding polynomial formalism to the const. pH mol. dynamics (CpHMD) framework, we are able to obtain the pH-dependent binding free energy profiles for various BACE-1-inhibitor complexes. Our results highlight the importance of correctly addressing the binding-induced protonation changes in protein-ligand systems where binding accompanies a net proton transfer. This work comprises the first application of our CpHMD-based free energy computational method to protein-ligand complexes and illustrates the value of CpHMD as an all-purpose tool for obtaining pH-dependent dynamics and binding free energies of biol. systems.
10
Sarkar, A.; Gupta, P. L.; Roitberg, A. E. pH-dependent conformational changes due to ionizable residues in a hydrophobic protein interior: The study of L25K and L125K variants of SNase. J. Phys. Chem. B 2019, 123, 5742– 5754, DOI: 10.1021/acs.jpcb.9b03816
Google Scholar
10
pH-Dependent Conformational Changes Due to Ionizable Residues in a Hydrophobic Protein Interior: The Study of L25K and L125K Variants of SNase
Sarkar, Ankita; Gupta, Pancham Lal; Roitberg, Adrian E.
Journal of Physical Chemistry B (2019), 123 (27), 5742-5754CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
Ionizable residues in the hydrophobic interior of certain proteins are known to play important roles in life processes like energy transduction and enzyme catalysis. These internal ionizable residues show exptl. apparent pKa values having large shifts as compared to their values in soln. In the present work, we study the pH-dependent conformational changes undergone by two variants of staphylococcal nuclease (SNase), L25K and L125K, using pH replica exchange mol. dynamics (pH-REMD) in explicit solvent. Our results show that the obsd. pKa of Lys25 and Lys125 are significantly different than their pKa in soln. We obsd. that the internal lysine residues prefer to be water-exposed when protonated at low pH, but they remain buried within the hydrophobic pocket when deprotonated at high pH. Using thermodn. laws, we est. the microscopic conformation-specific pKa of the water-exposed and buried conformations of the internal lysine residues and explain their relation to the macroscopic obsd. pKa values. We present the differences in the microscopic mechanisms that lead to similar exptl. obsd. apparent pKa of Lys25 and Lys125, and explain the need of thermodn. models of different complexities to account for our calcns. We see that L25K displays pH-dependent fluctuations throughout the entire β barrel and the α1 helix. In contrast, pH-independent fluctuations are obsd. in L125K, primarily limited to the α3 helix. The present computational study offers a detailed atomistic understanding of the determinants of the obsd. anomalous pKa of internal ionizable residues, bolstering the exptl. findings.
11
Sarkar, A.; Roitberg, A. E. PH-Dependent Conformational Changes Lead to a Highly Shifted p K a for a Buried Glutamic Acid Mutant of SNase. J. Phys. Chem. B 2020, 124, 11072– 11080, DOI: 10.1021/acs.jpcb.0c07136
Google Scholar
11
pH-Dependent Conformational Changes Lead to a Highly Shifted pKa for a Buried Glutamic Acid Mutant of SNase
Sarkar, Ankita; Roitberg, Adrian E.
Journal of Physical Chemistry B (2020), 124 (49), 11072-11080CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
Ionizable residues are rarely present in the hydrophobic interior of proteins, but when they are, they play important roles in biol. processes such as energy transduction and enzyme catalysis. Internal ionizable residues have anomalous exptl. pKa values with respect to their pKa in bulk water. This work studies the atomistic cause of the highly shifted pKa of the internal Glu23 in the artificially mutated variant V23E of Staphylococcal Nuclease (SNase) using pH replica exchange mol. dynamics (pH-REMD) simulations. The pKa of Glu23 obtained from the authors' calcns. is 6.55, which is elevated with respect to the glutamate pKa of 4.40 in bulk water. The calcd. value is close to the exptl. pKa of 7.10. The authors' simulations show that the highly shifted pKa of Glu23 is the product of a pH-dependent conformational change, which was obsd. exptl. and also seen in the authors' simulations. The authors carry out an anal. of this pH-dependent conformational change in response to the protonation state change of Glu23. Using a four-state thermodn. model, the authors est. the two conformation-specific pKa values of Glu23 and describe the coupling between the conformational and ionization equil.
12
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: Struct., Funct., Bioinf. 1997, 27, 523– 544, DOI: 10.1002/(SICI)1097-0134(199704)27:4<523::AID-PROT6>3.0.CO;2-B
Google Scholar
12
Simulation of protein conformational freedom as a function of pH: constant-pH molecular dynamics using implicit titration
Baptista, Antonio M.; Martel, Paulo J.; Petersen, Steffen B.
Proteins: Structure, Function, and Genetics (1997), 27 (4), 523-544CODEN: PSFGEY; ISSN:0887-3585. (Wiley-Liss)
Soln. pH is a determinant parameter on protein function and stability, and its inclusion in mol. dynamics simulations is attractive for studies at the mol. level. Current mol. dynamics simulations can consider pH only in a very limited way, through a somewhat arbitrary choice of a set of fixed charges on the titrable sites. Conversely, continuum electrostatic methods that explicitly treat pH effects assume a single protein conformation whose choice is not clearly defined. In this paper we describe a general method that combines both titrn. and conformational freedom. The method is based on a potential of mean force for implicit titrn. and combines both usual mol. dynamics and pH-dependent calcns. based on continuum methods. A simple implementation of the method, using a mean field approxn., is presented and applied to the bovine pancreatic trypsin inhibitor. We believe that this const.-pH mol. dynamics method, by correctly sampling both charges and conformation, can become a valuable help in the understanding of the dependence of protein function and stability on pH.
13
Baptista, A. M.; Teixeira, V. H.; Soares, C. M. Constant-pH molecular dynamics using stochastic titration. J. Chem. Phys. 2002, 117, 4184– 4200, DOI: 10.1063/1.1497164
Google Scholar
13
Constant-pH molecular dynamics using stochastic titration
Baptista, Antonio M.; Teixeira, Vitor H.; Soares, Claudio M.
Journal of Chemical Physics (2002), 117 (9), 4184-4200CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
A new method is proposed for performing const.-pH mol. dynamics (MD) simulations, i.e., MD simulations where pH is one of the external thermodn. parameters, like the temp. or the pressure. The protonation state of each titrable site in the solute is allowed to change during a mol. mechanics (MM) MD simulation, the new states being obtained from a combination of continuum electrostatics (CE) calcns. and Monte Carlo (MC) simulation of protonation equil. The coupling between the MM/MD and CE/MC algorithms is done in a way that ensures a proper Markov chain, sampling from the intended semigrand canonical distribution. This stochastic titrn. method is applied to succinic acid, aimed at illustrating the method and examg. the choice of its adjustable parameters. The complete titrn. of succinic acid, using const.-pH MD simulations at different pH values, gives a clear picture of the coupling between the trans/gauche isomerization and the protonation process, making it possible to reconcile some apparently contradictory results of previous studies. The present const.-pH MD method is shown to require a moderate increase of computational cost when compared to the usual MD method.
14
Bürgi, R.; Kollman, P. A.; van Gunsteren, W. F. Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation. Proteins: Struct., Funct., Bioinf. 2002, 47, 469– 480, DOI: 10.1002/prot.10046
Google Scholar
14
Simulating proteins at constant pH: an approach combining molecular dynamics and Monte Carlo simulation
Burgi, Roland; Kollman, Peter A.; Van Gunsteren, Wilfred F.
Proteins: Structure, Function, and Genetics (2002), 47 (4), 469-480CODEN: PSFGEY; ISSN:0887-3585. (Wiley-Liss, Inc.)
For the structure and function of proteins, the pH of the soln. is one of the detg. parameters. Current mol. dynamics (MD) simulations account for the soln. pH only in a limited way by keeping each titratable site in a chosen protonation state. We present an algorithm that generates trajectories at a Boltzmann distributed ensemble of protonation states by a combination of MD and Monte Carlo (MC) simulation. The algorithm is useful for pH-dependent structural studies and to investigate in detail the titrn. behavior of proteins. The method is tested on the acidic residues of the protein hen egg white lysozyme. It is shown that small structural changes may have a big effect on the pKA values of titratable residues.
15
Mongan, J.; Case, D. A.; McCammon, J. A. Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem. 2004, 25, 2038– 2048, DOI: 10.1002/jcc.20139
Google Scholar
15
Constant pH molecular dynamics in generalized Born implicit solvent
Mongan, John; Case, David A.; McCammon, J. Andrew
Journal of Computational Chemistry (2004), 25 (16), 2038-2048CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
A new method is proposed for const. pH mol. dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB-derived energies. The method is applied to four different crystal structures of hen egg-white lysozyme (HEWL). PKa predictions derived from the simulations have root-mean-square (RMS) error of 0.82 relative to exptl. values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics-only models. A strong correlation between conformation and protonation state is noted and quant. analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics.
16
Meng, Y.; Roitberg, A. E. Constant pH Replica Exchange Molecular Dynamics in Biomolecules Using a Discrete Protonation Model. J. Chem. Theory Comput. 2010, 6, 1401– 1412, DOI: 10.1021/ct900676b
Google Scholar
16
Constant pH Replica Exchange Molecular Dynamics in Biomolecules Using a Discrete Protonation Model
Meng, Yilin; Roitberg, Adrian E.
Journal of Chemical Theory and Computation (2010), 6 (4), 1401-1412CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
A const. pH replica exchange mol. dynamics (REMD) method is proposed and implemented to improve coupled protonation and conformational state sampling. By mixing conformational sampling at const. pH (with discrete protonation states) with a temp. ladder, this method avoids conformational trapping. Our method was tested and applied to seven different biol. systems. The const. pH REMD not only predicted pKa correctly for small, model compds. but also converged faster than const. pH mol. dynamics (MD). We further tested our const. pH REMD on a heptapeptide from the ovomucoid third domain (OMTKY3). Although const. pH REMD and MD produced very close pKa values, the const. pH REMD showed its advantage in the efficiency of conformational and protonation state samplings.
17
Itoh, S. G.; Damjanovic, A.; Brooks, B. R. pH replica-exchange method based on discrete protonation states. Proteins: Struct., Funct., Bioinf. 2011, 79, 3420– 3436, DOI: 10.1002/prot.23176
Google Scholar
17
pH replica-exchange method based on discrete protonation states
Itoh, Satoru G.; Damjanovic, Ana; Brooks, Bernard R.
Proteins: Structure, Function, and Bioinformatics (2011), 79 (12), 3420-3436CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
The authors propose a new algorithm for obtaining proton titrn. curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the const. pH algorithm of J. Mongan et al. (2004). In the original replica-exchange method, simulations of different replicas are performed at different temps., and the temps. are exchanged between the replicas. In the authors' PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original const. pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coeffs. of titratable residues. By performing multiple sets of const. pH simulations started with different initial states, the accuracy of predicted pKa values and Hill coeffs. obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titrn. points and thus better precision of measured pKa values and Hill coeffs. In addn., PHREM exhibits faster convergence of individual simulations than the original const. pH algorithm. Proteins 2011; © 2011 Wiley-Liss, Inc.
18
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, DOI: 10.1021/ct401042b
Google Scholar
18
Constant pH Replica Exchange Molecular Dynamics in Explicit Solvent Using Discrete Protonation States: Implementation, Testing, and Validation
Swails, Jason M.; York, Darrin M.; Roitberg, Adrian E.
Journal of Chemical Theory and Computation (2014), 10 (3), 1341-1352CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
By utilizing Graphics Processing Units, we show that const. pH mol. dynamics simulations (CpHMD) run in Generalized Born (GB) implicit solvent for long time scales can yield poor pKa predictions as a result of sampling unrealistic conformations. To address this shortcoming, we present a method for performing const. pH mol. dynamics simulations (CpHMD) in explicit solvent using a discrete protonation state model. The method involves std. mol. dynamics (MD) being propagated in explicit solvent followed by protonation state changes being attempted in GB implicit solvent at fixed intervals. Replica exchange along the pH-dimension (pH-REMD) helps to obtain acceptable titrn. behavior with the proposed method. We analyzed the effects of various parameters and settings on the titrn. behavior of CpHMD and pH-REMD in explicit solvent, including the size of the simulation unit cell and the length of the relaxation dynamics following protonation state changes. We tested the method with the amino acid model compds., a small pentapeptide with two titratable sites, and hen egg white lysozyme (HEWL). The proposed method yields superior predicted pKa values for HEWL over hundreds of nanoseconds of simulation relative to corresponding predicted values from simulations run in implicit solvent.
19
Swails, J. M.; Roitberg, A. E. Enhancing conformation and protonation state sampling of hen egg white lysozyme using pH replica exchange molecular dynamics. J. Chem. Theory Comput. 2012, 8, 4393– 4404, DOI: 10.1021/ct300512h
Google Scholar
19
Enhancing Conformation and Protonation State Sampling of Hen Egg White Lysozyme Using pH Replica Exchange Molecular Dynamics
Swails, Jason M.; Roitberg, Adrian E.
Journal of Chemical Theory and Computation (2012), 8 (11), 4393-4404CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
We evaluate the efficiency of the pH replica exchange mol. dynamics (pH-REMD) method proposed by Itoh et al. (Proteins 2011, 79, 3420-3436) by using it to predict the pKa values of the titratable residues in hen egg white lysozyme (HEWL). PKa values predicted using pH-REMD converge significantly faster than those calcd. using const. pH mol. dynamics (CpHMD). Furthermore, increasing the frequency between exchange attempts in pH-REMD simulations improves protonation and conformational state sampling. By enabling the simulation to sample both conformational and protonation states more rapidly, pH-REMD simulations provide valuable insight into the pH-dependence of HEWL that the CpHMD simulations failed to capture. We present an efficient and highly scalable implementation of pH-REMD as an attractive enhancement to traditional CpHMD methods.
20
Börjesson, U.; Hünenberger, P. H. Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines. J. Chem. Phys. 2001, 114, 9706– 9719, DOI: 10.1063/1.1370959
Google Scholar
20
Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines
Borjesson, Ulf; Hunenberger, Philippe H.
Journal of Chemical Physics (2001), 114 (22), 9706-9719CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
A method is developed for performing classical explicit-solvent mol. dynamics (MD) simulations at const. pH, where the protonation state of each ionizable (titratable) group in a simulated compd. is allowed to fluctuate in time, depending on the instantaneous system configuration and the imposed pH. In this method, each ionizable group is treated as a mixed state, i.e., the interaction-function parameters for the group are a linear combination of those of the protonated state and those of the deprotonated state. Free protons are not handled explicitly. Instead, the extent of deprotonation of each group is relaxed towards its equil. value by weak coupling to a "proton bath. " The method relies on precalibrated empirical functions, one for each type of ionizable group present in the simulated compd., which are obtained through multiple MD simulations of monofunctional model compds. In this study, the method is described in detail and its application illustrated by a series of const.-pH MD simulations of small monofunctional amines. In particular, we investigate the influence of the relaxation time used in the weak-coupling scheme, the choice of appropriate model compds. for the calibration of the required empirical functions, and corrections for finite-size effects linked with the small size of the simulation box.
21
Lee, M. S.; Salsbury, F. R., Jr.; Brooks, C. L., III Constant-pH molecular dynamics using continuous titration coordinates. Proteins: Struct., Funct., Bioinf. 2004, 56, 738– 752, DOI: 10.1002/prot.20128
Google Scholar
21
Constant-pH molecular dynamics using continuous titration coordinates
Lee, Michael S.; Salsbury, Freddie R., Jr.; Brooks, Charles L., III
Proteins: Structure, Function, and Bioinformatics (2004), 56 (4), 738-752CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
In this work, we explore the question of whether pKa calcns. based on a microscopic description of the protein and a macroscopic description of the solvent can be implemented to examine conformationally dependent proton shifts in proteins. To this end, we introduce a new method for performing const.-pH mol. dynamics (PHMD) simulations utilizing the generalized Born implicit solvent model. This approach employs an extended Hamiltonian in which continuous titrn. coordinates propagate simultaneously with the at. motions of the system. The values adopted by these coordinates are modulated by potentials of mean force of isolated titratable model groups and the pH to control the proton occupation at particular sites in the polypeptide. Our results for four different proteins yield an abs. av. error of ∼1.6 pK units, and point to the role that thermally driven relaxation of the protein environment in the vicinity of titrating groups plays in modulating the local pKa, thereby influencing the obsd. pK1/2 values. While the accuracy of our method is not yet equiv. to methods that obtain pK1/2 values through the ad hoc scaling of electrostatics, the present approach and const. pH methods in general provide a useful framework for studying pH-dependent phenomena. Further work to improve our model to approach quant. agreement with expt. is outlined.
22
Khandogin, J.; Brooks, C. L. Constant pH Molecular Dynamics with proton tautomerism. Biophys. J. 2005, 89, 141– 157, DOI: 10.1529/biophysj.105.061341
Google Scholar
22
Constant pH molecular dynamics with proton tautomerism
Khandogin, Jana; Brooks, Charles L., III
Biophysical Journal (2005), 89 (1), 141-157CODEN: BIOJAU; ISSN:0006-3495. (Biophysical Society)
The current article describes a new two-dimensional λ-dynamics method to include proton tautomerism in continuous const. pH mol. dynamics (CPHMD) simulations. The two-dimensional λ-dynamics framework is used to devise a tautomeric state titrn. model for the CPHMD simulations involving carboxyl and histidine residues. Combined with the GBSW implicit solvent model, the new method is tested on titrn. simulations of blocked histidine and aspartic acid as well as two benchmark proteins, turkey ovomucoid third domain (OMTKY3) and RNase A (RNase A). A detailed anal. of the errors inherent to the CPHMD methodol. as well as those due to the underlying solvation model is given. The av. abs. error for the computed pKa values in OMTKY3 is 1.0 pK unit. In RNase A the av. abs. errors for the carboxyl and histidine residues are 1.6 and 0.6 pK units, resp. In contrast to the previous work, the new model predicts the correct sign for all the pKa shifts, but one, in the benchmark proteins. The predictions of the tautomeric states of His12 and His48 and the conformational states of His48 and His119 are in agreement with expt. Based on the simulations of OMTKY3 and RNase A, the current work has demonstrated the capability of the CPHMD technique in revealing pH-coupled conformational dynamics of protein side chains.
23
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, DOI: 10.1021/ct200061r
Google Scholar
23
Constant pH Molecular Dynamics in Explicit Solvent with λ-Dynamics
Donnini, Serena; Tegeler, Florian; Groenhof, Gerrit; Grubmuller, Helmut
Journal of Chemical Theory and Computation (2011), 7 (6), 1962-1978CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
PH is an important parameter in condensed-phase systems, because it dets. the protonation state of titratable groups and thus influences the structure, dynamics, and function of mols. in soln. In most force field simulation protocols, however, the protonation state of a system (rather than its pH) is kept fixed and cannot adapt to changes of the local environment. Here, we present a method, implemented within the MD package GROMACS, for const. pH mol. dynamics simulations in explicit solvent that is based on the λ-dynamics approach. In the latter, the dynamics of the titrn. coordinate λ, which interpolates between the protonated and deprotonated states, is driven by generalized forces between the protonated and deprotonated states. The hydration free energy, as a function of pH, is included to facilitate const. pH simulations. The protonation states of titratable groups are allowed to change dynamically during a simulation, thus reproducing av. protonation probabilities at a certain pH. The accuracy of the method is tested against titrn. curves of single amino acids and a dipeptide in explicit solvent.
24
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, DOI: 10.1021/ct200146j
Google Scholar
24
Continuous Constant pH Molecular Dynamics in Explicit Solvent with pH-Based Replica Exchange
Wallace, Jason A.; Shen, Jana K.
Journal of Chemical Theory and Computation (2011), 7 (8), 2617-2629CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
A computational tool that offers accurate pKa values and atomically detailed knowledge of protonation-coupled conformational dynamics is valuable for elucidating mechanisms of energy transduction processes in biol., such as enzyme catalysis and electron transfer as well as proton and drug transport. Toward this goal the authors present a new technique of embedding continuous const. pH mol. dynamics within an explicit-solvent representation. In this technique the authors make use of the efficiency of the generalized-Born (GB) implicit-solvent model for estg. the free energy of protein solvation while propagating conformational dynamics using the more accurate explicit-solvent model. Also, the authors employ a pH-based replica exchange scheme to significantly enhance both protonation and conformational state sampling. Benchmark data of five proteins including HP36, NTL9, BBL, HEWL, and SNase yield an av. abs. deviation of 0.53 and a root mean squared deviation of 0.74 from exptl. data. This level of accuracy is obtained with 1 ns simulations per replica. Detailed anal. reveals that explicit-solvent sampling provides increased accuracy relative to the previous GB-based method by preserving the native structure, providing a more realistic description of conformational flexibility of the hydrophobic cluster, and correctly modeling solvent mediated ion-pair interactions. Thus, the authors anticipate that the new technique will emerge as a practical tool to capture ionization equil. while enabling an intimate view of ionization coupled conformational dynamics that is difficult to delineate with exptl. techniques alone.
25
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, DOI: 10.1063/1.4766352
Google Scholar
25
Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH
Wallace, Jason A.; Shen, Jana K.
Journal of Chemical Physics (2012), 137 (18), 184105/1-184105/9CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
Recent development of const. pH mol. dynamics (CpHMD) methods has offered promise for adding pH-stat in mol. dynamics simulations. However, until now the working pH mol. dynamics (pHMD) implementations are dependent in part or whole on implicit-solvent models. Here we show that proper treatment of long-range electrostatics and maintaining charge neutrality of the system are crit. for extending the continuous pHMD framework to the all-atom representation. The former is achieved here by adding forces to titrn. coordinates due to long-range electrostatics based on the generalized reaction field method, while the latter is made possible by a charge-leveling technique that couples proton titrn. with simultaneous ionization or neutralization of a co-ion in soln. We test the new method using the pH-replica-exchange CpHMD simulations of a series of aliph. dicarboxylic acids with varying carbon chain length. The av. abs. deviation from the exptl. pKa values is merely 0.18 units. The results show that accounting for the forces due to extended electrostatics removes the large random noise in propagating titrn. coordinates, while maintaining charge neutrality of the system improves the accuracy in the calcd. electrostatic interaction between ionizable sites. Thus, we believe that the way is paved for realizing pH-controlled all-atom mol. dynamics in the near future. (c) 2012 American Institute of Physics.
26
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, DOI: 10.1021/ct2006314
Google Scholar
26
Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent
Goh, Garrett B.; Knight, Jennifer L.; Brooks, Charles L.
Journal of Chemical Theory and Computation (2012), 8 (1), 36-46CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
The nucleosides of adenine and cytosine have pKa values of 3.50 and 4.08, resp., and are assumed to be unprotonated under physiol. conditions. However, evidence from recent NMR and X-ray crystallog. studies has revealed the prevalence of protonated adenine and cytosine in RNA macromols. Such nucleotides with elevated pKa values may play a role in stabilizing RNA structure and participate in the mechanism of ribozyme catalysis. With the work presented here, we establish the framework and demonstrate the first const. pH MD simulations (CPHMD) for nucleic acids in explicit solvent, in which the protonation state is coupled to the dynamical evolution of the RNA system via λ-dynamics. We adopt the new functional form λNexp for λ that was recently developed for multisite λ-dynamics (MSλD) and demonstrate good sampling characteristics in which rapid and frequent transitions between the protonated and unprotonated states at pH = pKa are achieved. Our calcd. pKa values of simple nucleotides are in a good agreement with exptl. measured values, with a mean abs. error of 0.24 pKa units. This work demonstrates that CPHMD can be used as a powerful tool to investigate pH-dependent biol. properties of RNA macromols.
27
Chen, W.; Wallace, J. A.; Yue, Z.; Shen, J. K. Introducing Titratable Water to All-Atom Molecular Dynamics at Constant pH. Biophys. J. 2013, 105, L15– L17, DOI: 10.1016/j.bpj.2013.06.036
Google Scholar
27
Introducing Titratable Water to All-Atom Molecular Dynamics at Constant pH
Chen, Wei; Wallace, Jason A.; Yue, Zhi; Shen, Jana K.
Biophysical Journal (2013), 105 (4), L15-L17CODEN: BIOJAU; ISSN:0006-3495. (Cell Press)
Recent development of titratable coions has paved the way for realizing all-atom mol. dynamics at const. pH. To further improve phys. realism, here we describe a technique in which proton titrn. of the solute is directly coupled to the interconversion between water and hydroxide or hydronium. We test the new method in replica-exchange continuous const. pH mol. dynamics simulations of three proteins, HP36, BBL, and HEWL. The calcd. pKa values based on 10-ns sampling per replica have the av. abs. and root-mean-square errors of 0.7 and 0.9 pH units, resp. Introducing titratable water in mol. dynamics offers a means to model proton exchange between solute and solvent, thus opening a door to gaining new insights into the intricate details of biol. phenomena involving proton translocation.
28
Lee, J.; Miller, B. T.; Damjanovic, A.; Brooks, B. R. Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange. J. Chem. Theory Comput. 2014, 10, 2738– 2750, DOI: 10.1021/ct500175m
Google Scholar
28
Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange
Lee, Juyong; Miller, Benjamin T.; Damjanovic, Ana; Brooks, Bernard R.
Journal of Chemical Theory and Computation (2014), 10 (7), 2738-2750CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
We present a new computational approach for const. pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike const. pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The anal. is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the min. energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic acid, lysine, and glutamic acid, which have two protonation states, a histidine with three states, a four-residue peptide with four states, and snake cardiotoxin with eight states. The pKa values estd. with the EDS-HREX method agree well with the exptl. pKa values. The mean abs. errors of small benchmark systems range from 0.03 to 0.17 pKa units, and those of three titratable groups of snake cardiotoxin range from 0.2 to 1.6 pKa units. This study demonstrates that EDS-HREX is a potent theor. framework, which gives the correct description of multiple protonation states and good calcd. pKa values.
29
Dobrev, P.; Donnini, S.; Groenhof, G.; Grubmüller, H. Accurate Three States Model for Amino Acids with Two Chemically Coupled Titrating Sites in Explicit Solvent Atomistic Constant pH Simulations and p K a Calculations. J. Chem. Theory Comput. 2017, 13, 147– 160, DOI: 10.1021/acs.jctc.6b00807
Google Scholar
29
Accurate Three States Model for Amino Acids with Two Chemically Coupled Titrating Sites in Explicit Solvent Atomistic Constant pH Simulations and pKa Calculations
Dobrev, Plamen; Donnini, Serena; Groenhof, Gerrit; Grubmueller, Helmut
Journal of Chemical Theory and Computation (2017), 13 (1), 147-160CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
Correct protonation of titratable groups in biomols. is crucial for their accurate description by mol. dynamics simulations. In the context of const. pH simulations, an addnl. protonation degree of freedom is introduced for each titratable site, allowing the protonation state to change dynamically with changing structure or electrostatics. Here, we implement a second reaction coordinate which switches between two tautomeric states of an amino acid with chem. coupled titratable sites, such as aspartate (Asp), glutamate (Glu), and histidine (His). To this aim, we test a scheme involving three protonation states. To facilitate charge neutrality as required for periodic boundary conditions and Particle Mesh Ewald (PME) electrostatics, titrn. of each resp. amino acid is coupled to a "water" mol. that is charged in opposite direction. Addnl., a force field modification for Amber99sb is introduced and tested for the description of carboxyl group protonation. Our three states model is tested by titrn. simulations of Asp, Glu, and His, yielding a good agreement, reproducing the correct geometry of the groups in their different protonation forms. We further show that the ion concn. change due to the neutralizing "water" mols. does not significantly affect the protonation free energies of the titratable groups, suggesting that the three states model provides a good description of biomol. dynamics at const. pH.
30
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, DOI: 10.1021/acs.jctc.6b00552
Google Scholar
30
All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water
Huang, Yandong; Chen, Wei; Wallace, Jason A.; Shen, Jana
Journal of Chemical Theory and Computation (2016), 12 (11), 5411-5421CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
Development of a pH stat to properly control soln. pH in biomol. simulations has been a long-standing goal in the community. Towards this goal recent years have witnessed the emergence of the so-called const. pH mol. dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or electrostatic truncation schemes. Here the authors report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous const. pH mol. dynamics (CpHMD) method, enabling CpHMD to be performed with a std. MD engine at a fractional added computational cost. The authors demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP38, BBL, HEWL and SNase. With the sampling time of 10 ns per replica, most pKa's are converged, yielding the av. abs. and root-mean-square deviations of 0.61 and 0.77, resp., from expt. Linear regression of the calcd. vs. exptl. pKa shifts gives a correlation coeff. of 0.79, a slope of 1 and an intercept near 0. Anal. reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromol. systems in various environments, enabling at. insights into pH-dependent phenomena involving not only sol. proteins but also transmembrane proteins, nucleic acids, surfactants and polysaccharides.
31
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, DOI: 10.1021/acs.jcim.8b00227
Google Scholar
31
Generalized Born Based Continuous Constant pH Molecular Dynamics in Amber: Implementation, Benchmarking and Analysis
Huang, Yandong; Harris, Robert C.; Shen, Jana
Journal of Chemical Information and Modeling (2018), 58 (7), 1372-1383CODEN: JCISD8; ISSN:1549-9596. (American Chemical Society)
Soln. pH plays an important role in structure and dynamics of biomol. systems; however, pH effects cannot be accurately accounted for in conventional mol. dynamics simulations based on fixed protonation states. Continuous const. pH mol. dynamics (CpHMD) based on the λ-dynamics framework calcs. protonation states on the fly during dynamical simulation at a specified pH condition. Here the authors report the CPU-based implementation of the CpHMD method based on the GBNeck2 generalized Born (GB) implicit-solvent model in the pmemd engine of the Amber mol. dynamics package. The performance of the method was tested using pH replica-exchange titrn. simulations of Asp, Glu and His side chains in 4 miniproteins and 7 enzymes with exptl. known pKa's, some of which are significantly shifted from the model values. The added computational cost due to CpHMD titrn. ranges from 11 to 33% for the data set and scales roughly linearly as the ratio between the titrable sites and no. of solute atoms. Comparison of the exptl. and calcd. pKa's using 2 ns per replica sampling yielded a mean unsigned error of 0.70, a root-mean-squared error of 0.91, and a linear correlation coeff. of 0.79. Though this level of accuracy is similar to the GBSW-based CpHMD in CHARMM, in contrast to the latter, the current implementation was able to reproduce the exptl. orders of the pKa's of the coupled carboxylic dyads. The authors quantified the sampling errors, which revealed that prolonged simulation is needed to converge pKa's of several titratable groups involved in salt-bridge-like interactions or deeply buried in the protein interior. The authors' benchmark data demonstrate that GBNeck2-CpHMD is an attractive tool for protein pKa predictions.
32
Harris, R. C.; Shen, J. GPU-accelerated implementation of continuous constant pH molecular dynamics in amber: pKa predictions with single-pH simulations. J. Chem. Inf. Model. 2019, 59, 4821– 4832, DOI: 10.1021/acs.jcim.9b00754
Google Scholar
32
GPU-Accelerated Implementation of Continuous Constant pH Molecular Dynamics in Amber: pKa Predictions with Single-pH Simulations
Harris, Robert C.; Shen, Jana
Journal of Chemical Information and Modeling (2019), 59 (11), 4821-4832CODEN: JCISD8; ISSN:1549-9596. (American Chemical Society)
The authors present a GPU implementation of the continuous const. pH mol. dynamics (CpHMD) based on the most recent generalized Born implicit-solvent model in the pmemd engine of the Amber mol. dynamics package. To test the accuracy of the tool for rapid pKa predictions, a series of 2-ns single-pH simulations were performed for over 120 titratable residues in 10 benchmark proteins that were previously used to test the various continuous CpHMD methods. The calcd. pKa's showed a root-mean-square deviation of 0.80 and correlation coeff. of 0.83 with respect to expt. 90% of the pKa's were converged with estd. errors <0.1 pH units. Surprisingly, this level of accuracy is similar to the authors' previous replica-exchange simulations with 2 ns per replica and an exchange attempt frequency of 2 ps-1 (Huang, Harris and Shen, J Chem Info Model, 2018). For the linked titrn. sites in two enzymes, although residue-specific protonation state sampling in the single-pH simulations was not converged within 2 ns, the protonation fraction of the linked residues appeared to be largely converged, and the exptl. macroscopic {\pka} values were reproduced to within 1 pH unit. Comparison with replica-exchange simulations with different exchange attempt frequencies showed that the splitting between the two macroscopic pKa's is underestimated with frequent exchange attempts such as 2 ps$̂{-1}$, while single-pH simulations overestimate the splitting. The same trend is seen for the single-pH vs. replica-exchange simulations of a hydrogen-bonded aspartyl dyad in a much larger protein. A 2-ns single-pH simulation of a 400-residue protein takes about one hour on a single NVIDIA GeForce RTX 2080 graphics card, which is over 1000 times faster than a CpHMD run on a single CPU core of a high-performance computing cluster node. Thus, the authors envision that GPU-accelerated continuous CpHMD may be used in routine pKa predictions for a variety of applications, from assisting MD simulations with protonation state assignment to offering pH-dependent corrections of binding free energies and identifying reactive hot spots for covalent drug design.
33
Harris, R. C.; Liu, R.; Shen, J. Predicting reactive cysteines with implicit-solvent-based continuous constant pH molecular dynamics in amber. J. Chem. Theory Comput. 2020, 16, 3689– 3698, DOI: 10.1021/acs.jctc.0c00258
Google Scholar
33
Predicting Reactive Cysteines with Implicit-Solvent-Based Continuous Constant pH Molecular Dynamics in Amber
Harris, Robert C.; Liu, Ruibin; Shen, Jana
Journal of Chemical Theory and Computation (2020), 16 (6), 3689-3698CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
Cysteines existing in the deprotonated thiolate form or having a tendency to become deprotonated are important players in enzymic and cellular redox functions and frequently exploited in covalent drug design; however, most computational studies assume cysteines as protonated. Thus, developing an efficient tool that can make accurate and reliable predictions of cysteine protonation states is timely needed. The authors recently implemented a generalized Born (GB) based continuous const. pH mol. dynamics (CpHMD) method in Amber for protein pKa calcns. on CPUs and GPUs. Here the authors benchmark the performance of GB-CpHMD for predictions of cysteine pKa's and reactivities using a data set of 24 proteins with both down- and upshifted cysteine pKa's. The authors found that 10 ns single-pH or 4 ns replica-exchange CpHMD titrns. gave root-mean-square errors of 1.2-1.3 and correlation coeffs. of 0.8-0.9 with respect to expt. The accuracy of predicting thiolates or reactive cysteines at physiol. pH with single-pH titrns. is 86 or 81% with a precision of 100 or 90%, resp. This performance well surpasses the traditional structure-based methods, particularly a widely used empirical pKa tool that gives an accuracy less than 50%. The authors discuss simulation convergence, dependence on starting structures, common determinants of the pKa downshifts and upshifts, and the origin of the discrepancies from the structure-based calcns. The work suggests that CpHMD titrns. can be performed on a desktop computer equipped with a single GPU card to predict cysteine protonation states for a variety of applications, from understanding biol. functions to covalent drug design.
34
Grünewald, F.; Souza, P. C.; Abdizadeh, H.; Barnoud, J.; de Vries, A. H.; Marrink, S. J. Titratable Martini model for constant pH simulations. J. Chem. Phys. 2020, 153, 024118, DOI: 10.1063/5.0014258
Google Scholar
34
Titratable Martini model for constant pH simulations
Gruenewald, Fabian; Souza, Paulo C. T.; Abdizadeh, Haleh; Barnoud, Jonathan; de Vries, Alex H.; Marrink, Siewert J.
Journal of Chemical Physics (2020), 153 (2), 024118CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
The authors deliver a proof of concept for a fast method that introduces pH effects into classical coarse-grained (CG) mol. dynamics simulations. The authors' approach is based upon the latest version of the popular Martini CG model to which explicit proton mimicking particles were added. The authors verify the authors' approach against exptl. data involving several different mols. and different environmental conditions. In particular, the authors compute titrn. curves, pH dependent free energies of transfer, and lipid bilayer membrane affinities as a function of pH. Using oleic acid as an example compd., the authors further illustrate that the authors' method can be used to study passive translocation in lipid bilayers via protonation. Finally, the authors' model reproduces qual. the expansion of the macromol. dendrimer poly(propylene imine) as well as the assocd. pKa shift of its different generations. This example demonstrates that the authors' model is able to pick up collective interactions between titratable sites in large mols. comprising many titratable functional groups. (c) 2020 American Institute of Physics.
35
Mongan, J.; Case, D. A. Biomolecular simulations at constant pH. Curr. Opin. Struct. Biol. 2005, 15, 157– 163, DOI: 10.1016/j.sbi.2005.02.002
Google Scholar
35
Biomolecular simulations at constant pH
Mongan, John; Case, David A.
Current Opinion in Structural Biology (2005), 15 (2), 157-163CODEN: COSBEF; ISSN:0959-440X. (Elsevier Ltd.)
A review. Like temp. and pressure, the soln. pH is an important thermodn. variable that is commonly varied in expts. and is used by cells to influence biochem. function. It is now becoming feasible to carry out practical mol. dynamics simulations that mimic the thermodn. of such expts., by allowing proton transfer between the system of interest and a hypothetical bath of protons at a given pH. These are demanding calcns., because the energetics of charge changes upon protonation or deprotonation must be accurately modeled, and because such simulations must sample both mol. configurations and the large no. of protonation states that are possible for a mol. with many titrating sites.
36
Damjanovic, A.; Miller, B. T.; Okur, A.; Brooks, B. R. Reservoir pH replica exchange. J. Chem. Phys. 2018, 149, 072321, DOI: 10.1063/1.5027413
Google Scholar
36
Reservoir pH replica exchange
Damjanovic, Ana; Miller, Benjamin T.; Okur, Asim; Brooks, Bernard R.
Journal of Chemical Physics (2018), 149 (7), 072321/1-072321/12CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
We present the reservoir pH replica exchange (R-pH-REM) method for const. pH simulations. The R-pH-REM method consists of a two-step procedure; the first step involves generation of one or more reservoirs of conformations. Each reservoir is obtained from a std. or enhanced mol. dynamics simulation with a constrained (fixed) protonation state. In the second step, fixed charge constraints are relaxed, as the structures from one or more reservoirs are periodically injected into a const. pH or a pH-replica exchange (pH-REM) simulation. The benefit of this two-step process is that the computationally intensive part of conformational search can be decoupled from const. pH simulations, and various techniques for enhanced conformational sampling can be applied without the need to integrate such techniques into the pH-REM framework. Simulations on blocked Lys, KK, and KAAE peptides were used to demonstrate an agreement between pH-REM and R-pH-REM simulations. While the reservoir simulations are not needed for these small test systems, the real need arises in cases when ionizable mols. can sample two or more conformations sepd. by a large energy barrier, such that adequate sampling is not achieved on a time scale of std. const. pH simulations. Such problems might be encountered in protein systems that exploit conformational transitions for function. A hypothetical case is studied, a small mol. with a large torsional barrier; while results of pH-REM simulations depend on the starting structure, R-pH-REM calcns. on this model system are in excellent agreement with a theor. model. (c) 2018 American Institute of Physics.
37
Chen, Y.; Roux, B. Constant-pH Hybrid Nonequilibrium Molecular Dynamics - Monte Carlo Simulation Method. J. Chem. Theory Comput. 2015, 11, 3919– 3931, DOI: 10.1021/acs.jctc.5b00261
Google Scholar
37
Constant-pH Hybrid Nonequilibrium Molecular Dynamics-Monte Carlo Simulation Method
Chen, Yunjie; Roux, Benoit
Journal of Chemical Theory and Computation (2015), 11 (8), 3919-3931CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
A computational method is developed to carry out explicit solvent simulations of complex mol. systems under conditions of const. pH. In const.-pH simulations, preidentified ionizable sites are allowed to spontaneously protonate and deprotonate as a function of time in response to the environment and the imposed pH. The method consists of carrying out short nonequil. mol. dynamics (neMD) switching trajectories to generate phys. plausible configurations with changed protonation states that are subsequently accepted or rejected according to a Metropolis Monte Carlo (MC) criterion. To ensure microscopic detailed balance arising from such nonequil. switches, the at. momenta are altered according to the sym. two-ends momentum reversal prescription. To achieve higher efficiency, the original neMD-MC scheme is sepd. into two steps, reducing the need for generating a large no. of unproductive and costly nonequil. trajectories. In the first step, the protonation state of a site is randomly attributed via a Metropolis MC process on the basis of an intrinsic pKa; an attempted nonequil. switch is generated only if this change in protonation state is accepted. This hybrid two-step inherent pKa neMD-MC simulation method is tested with single amino acids in soln. (Asp, Glu, and His) and then applied to turkey ovomucoid third domain and hen egg-white lysozyme. Because of the simple linear increase in the computational cost relative to the no. of titratable sites, the present method is naturally able to treat extremely large systems.
38
Kong, X.; Brooks, C. L., III λ-dynamics: A new approach to free energy calculations. J. Chem. Phys. 1996, 105, 2414– 2423, DOI: 10.1063/1.472109
Google Scholar
There is no corresponding record for this reference.
39
Simonson, T.; Carlsson, J.; Case, D. A. Proton binding to proteins: p K a calculations with explicit and implicit solvent models. J. Am. Chem. Soc. 2004, 126, 4167– 4180, DOI: 10.1021/ja039788m
Google Scholar
39
Proton Binding to Proteins: pKa Calculations with Explicit and Implicit Solvent Models
Simonson, Thomas; Carlsson, Jens; Case, David A.
Journal of the American Chemical Society (2004), 126 (13), 4167-4180CODEN: JACSAT; ISSN:0002-7863. (American Chemical Society)
Ionizable residues play important roles in protein structure and activity, and proton binding is a valuable reporter of electrostatic interactions in these systems. The authors use mol. dynamics free energy simulations (MDFE) to compute proton pKa shifts, relative to a model compd. in soln., for three aspartate side chains in two proteins. Simulations with explicit solvent and with an implicit, dielec. continuum solvent are reported. The implicit solvent simulations use the generalized Born (GB) model, which provides an approx., anal. soln. to Poisson's equation. With explicit solvent, the direction of the pKa shifts is correct in all three cases with one force field (AMBER) and in two out of three cases with another (CHARMM). For two aspartates, the dielec. response to ionization is found to be linear, even though the sep. protein and solvent responses can be nonlinear. For thioredoxin Asp-26, nonlinearity arises from the presence of two substates that correspond to the two possible orientations of the protonated carboxylate. For this side chain, which is partly buried and has a large pKa upshift, very long simulations are needed to correctly sample several slow degrees of freedom that reorganize in response to the ionization. Thus, nearby Lys-57 rotates to form a salt bridge and becomes buried, while three waters intercalate along the opposite edge of Asp-26. Such strong and anisotropic reorganization is very difficult to predict with Poisson-Boltzmann methods that only consider electrostatic interactions and employ a single protein structure. In contrast, MDFE with a GB dielec. continuum solvent, used for the first time for pKa calcns., can describe protein reorganization accurately and gives encouraging agreement with expt. and with the explicit solvent simulations.
40
Chen, W.; Shen, J. K. Effects of System Net Charge and Electrostic Truncation on All-Atom Constant pH Molecular Dynamics. J. Comput. Chem. 2014, 35, 1986– 1996, DOI: 10.1002/jcc.23713
Google Scholar
40
Effects of system net charge and electrostatic truncation on all-atom constant pH molecular dynamics
Chen, Wei; Shen, Jana K.
Journal of Computational Chemistry (2014), 35 (27), 1986-1996CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
Const. pH mol. dynamics offers a means to rigorously study the effects of soln. pH on dynamical processes. Here, we address two crit. questions arising from the most recent developments of the all-atom continuous const. pH mol. dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states. (2) Is the enforcement of elec. neutrality necessary for const. pH simulations. We first examd. how the generalized reaction field and force-shifting schemes modify the electrostatic forces on the titrn. coordinates. Free energy simulations of model compds. were then carried out to delineate the errors in the deprotonation free energy and salt-bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titrn. of a mini-protein HP36 was used to understand the manifestation of the two types of errors in the calcd. pKa values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via co-titrating ions significantly improves the accuracy of protonation-state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge-compensating background plasma.
41
Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812– 7824, DOI: 10.1021/jp071097f
Google Scholar
41
The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations
Marrink, Siewert J.; Risselada, H. Jelger; Yefimov, Serge; Tieleman, D. Peter; De Vries, Alex H.
Journal of Physical Chemistry B (2007), 111 (27), 7812-7824CODEN: JPCBFK; ISSN:1520-6106. (American Chemical Society)
We present an improved and extended version of our coarse grained lipid model. The new version, coined the MARTINI force field, is parametrized in a systematic way, based on the reprodn. of partitioning free energies between polar and apolar phases of a large no. of chem. compds. To reproduce the free energies of these chem. building blocks, the no. of possible interaction levels of the coarse-grained sites has increased compared to those of the previous model. Application of the new model to lipid bilayers shows an improved behavior in terms of the stress profile across the bilayer and the tendency to form pores. An extension of the force field now also allows the simulation of planar (ring) compds., including sterols. Application to a bilayer/cholesterol system at various concns. shows the typical cholesterol condensation effect similar to that obsd. in all atom representations.
42
Donnini, S.; Ullmann, R. T.; Groenhof, G.; Grubmüller, H. Charge-neutral constant pH molecular dynamics simulations using a parsimonious proton buffer. J. Chem. Theory Comput. 2016, 12, 1040– 1051, DOI: 10.1021/acs.jctc.5b01160
Google Scholar
42
Charge-Neutral Constant pH Molecular Dynamics Simulations Using a Parsimonious Proton Buffer
Donnini, Serena; Ullmann, R. Thomas; Groenhof, Gerrit; Grubmuller, Helmut
Journal of Chemical Theory and Computation (2016), 12 (3), 1040-1051CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
In const. pH mol. dynamics simulations, the protonation states of titratable sites can respond to changes of the pH and of their electrostatic environment. Consequently, the no. of protons bound to the biomol., and therefore the overall charge of the system, fluctuates during the simulation. To avoid artifacts assocd. with a non-neutral simulation system, we introduce an approach to maintain neutrality of the simulation box in const. pH mol. dynamics simulations, while maintaining an accurate description of all protonation fluctuations. Specifically, we introduce a proton buffer that, like a buffer in expt., can exchange protons with the biomol. enabling its charge to fluctuate. To keep the total charge of the system const., the uptake and release of protons by the buffer are coupled to the titrn. of the biomol. with a constraint. We find that, because the fluctuation of the total charge (no. of protons) of a typical biomol. is much smaller than the no. of titratable sites of the biomol., the no. of buffer sites required to maintain overall charge neutrality without compromising the charge fluctuations of the biomol., is typically much smaller than the no. of titratable sites, implying markedly enhanced simulation and sampling efficiency.
43
Dobrev, P.; Vemulapalli, S. P. B.; Nath, N.; Griesinger, C.; Grubmüller, H. Probing the accuracy of explicit solvent constant pH molecular dynamics simulations for peptides. J. Chem. Theory Comput. 2020, 16, 2561– 2569, DOI: 10.1021/acs.jctc.9b01232
Google Scholar
43
Probing the Accuracy of Explicit Solvent Constant pH Molecular Dynamics Simulations for Peptides
Dobrev, Plamen; Vemulapalli, Sahithya Phani Babu; Nath, Nilamoni; Griesinger, Christian; Grubmueller, Helmut
Journal of Chemical Theory and Computation (2020), 16 (4), 2561-2569CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
Protonation states of titratable amino acids play a key role in many biomol. processes. Knowledge of protonatable residue charges at a given pH is essential for a correct understanding of protein catalysis, inter- and intramol. interactions, substrate binding, and protein dynamics for instance. However, acquiring exptl. values for individual amino acid protonation states of complex systems is not straightforward; therefore, several in silico approaches have been developed to tackle this issue. In this work, the authors assess the accuracy of the previously developed const. pH MD approach by comparing the theor. obtained pKa values for titratable residues with exptl. values from an equiv. NMR study. The authors selected a set of four pentapeptides, of adequately small size to ensure comprehensive sampling, but concurrently, due to their charge compn., posing a challenge for protonation state calcn. The comparison of the pKa values shows good agreement of the exptl. and the theor. approach with a largest difference of 0.25 pKa units. Further, the corresponding titrn. curves are in fair agreement, although the shift of the Hill coeff. from a value of 1 was not always reproduced in simulations. The phase space overlap in Cartesian space between trajectories generated in const. pH and std. MD simulations is fair and suggests that the const. pH MD approach reasonably well preserves the dynamics of the system, allowing dynamic protonation MD simulations without introducing structural artifacts.
44
Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089– 10092, DOI: 10.1063/1.464397
Google Scholar
44
Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems
Darden, Tom; York, Darrin; Pedersen, Lee
Journal of Chemical Physics (1993), 98 (12), 10089-92CODEN: JCPSA6; ISSN:0021-9606.
An N·log(N) method for evaluating electrostatic energies and forces of large periodic systems is presented. The method is based on interpolation of the reciprocal space Ewald sums and evaluation of the resulting convolution using fast Fourier transforms. Timings and accuracies are presented for three large cryst. ionic systems.
45
Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577– 8593, DOI: 10.1063/1.470117
Google Scholar
45
A smooth particle mesh Ewald method
Essmann, Ulrich; Perera, Lalith; Berkowitz, Max L.; Darden, Tom; Lee, Hsing; Pedersen, Lee G.
Journal of Chemical Physics (1995), 103 (19), 8577-93CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
The previously developed particle mesh Ewald method is reformulated in terms of efficient B-spline interpolation of the structure factors. This reformulation allows a natural extension of the method to potentials of the form 1/rp with p ≥ 1. Furthermore, efficient calcn. of the virial tensor follows. Use of B-splines in the place of Lagrange interpolation leads to analytic gradients as well as a significant improvement in the accuracy. The authors demonstrate that arbitrary accuracy can be achieved, independent of system size N, at a cost that scales as N log(N). For biomol. systems with many thousands of atoms and this method permits the use of Ewald summation at a computational cost comparable to that of a simple truncation method of 10 Å or less.
46
Knight, J. L.; Brooks, C. L., III Multisite λ dynamics for simulated structure-activity relationship studies. J. Chem. Theory Comput. 2011, 7, 2728– 2739, DOI: 10.1021/ct200444f
Google Scholar
46
Multisite λ Dynamics for Simulated Structure-Activity Relationship Studies
Knight, Jennifer L.; Brooks, Charles L., III
Journal of Chemical Theory and Computation (2011), 7 (9), 2728-2739CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
Multisite λ dynamics (MSλD) is a new free energy simulation method that is based on λ dynamics. It has been developed to enable multiple substituents at multiple sites on a common ligand core to be modeled simultaneously and their free energies assessed. The efficacy of MSλD for estg. relative hydration free energies and relative binding affinties is demonstrated using three test systems. Model compds. representing multiple identical benzene, dihydroxybenzene, and dimethoxybenzene mols. show that total combined MSλD trajectory lengths of ∼1.5 ns are sufficient to reliably achieve relative hydration free energy ests. within 0.2 kcal/mol and are less sensitive to the no. of trajectories that are used to generate these ests. for hybrid ligands that contain up to 10 substituents modeled at a single site or five substituents modeled at each of two sites. Relative hydration free energies among six benzene derivs. calcd. from MSλD simulations are in very good agreement with those from alchem. free energy simulations (with av. unsigned differences of 0.23 kcal/mol and R2 = 0.991) and the expt. (with av. unsigned errors of 1.8 kcal/mol and R2 = 0.959). Ests. of the relative binding affinities among 14 inhibitors of HIV-1 reverse transcriptase obtained from MSλD simulations are in reasonable agreement with those from traditional free energy simulations and the expt. (av. unsigned errors of 0.9 kcal/mol and R2 = 0.402). For the same level of accuracy and precision, MSλD simulations are achieved ∼20-50 times faster than traditional free energy simulations and thus with reliable force field parameters can be used effectively to screen tens to hundreds of compds. in structure-based drug design applications.
47
Goh, G. B.; Hulbert, B. S.; Zhou, H.; Brooks, C. L., III Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism. Proteins: Struct., Funct., Bioinf. 2014, 82, 1319– 1331, DOI: 10.1002/prot.24499
Google Scholar
47
Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism
Goh, Garrett B.; Hulbert, Benjamin S.; Zhou, Huiqing; Brooks, Charles L., III
Proteins: Structure, Function, and Bioinformatics (2014), 82 (7), 1319-1331CODEN: PSFBAF ISSN:. (Wiley-Blackwell)
PH is a ubiquitous regulator of biol. activity, including protein-folding, protein-protein interactions, and enzymic activity. Existing const. pH mol. dynamics (CPHMD) models that were developed to address questions related to the pH-dependent properties of proteins are largely based on implicit solvent models. However, implicit solvent models are known to underestimate the desolvation energy of buried charged residues, increasing the error assocd. with predictions that involve internal ionizable residue that are important in processes like hydrogen transport and electron transfer. Furthermore, discrete water and ions cannot be modeled in implicit solvent, which are important in systems like membrane proteins and ion channels. We report on an explicit solvent const. pH mol. dynamics framework based on multi-site λ-dynamics (CPHMDMSλD). In the CPHMDMSλD framework, we performed seamless alchem. transitions between protonation and tautomeric states using multi-site λ-dynamics, and designed novel biasing potentials to ensure that the phys. end-states are predominantly sampled. We show that explicit solvent CPHMDMSλD simulations model realistic pH-dependent properties of proteins such as the Hen-Egg White Lysozyme (HEWL), binding domain of 2-oxoglutarate dehydrogenase (BBL) and N-terminal domain of ribosomal protein L9 (NTL9), and the pKa predictions are in excellent agreement with exptl. values, with a RMSE ranging from 0.72 to 0.84 pKa units. With the recent development of the explicit solvent CPHMDMSλD framework for nucleic acids, accurate modeling of pH-dependent properties of both major class of biomols.-proteins and nucleic acids is now possible.
48
Knight, J. L.; Brooks, C. L., III Applying efficient implicit nongeometric constraints in alchemical free energy simulations. Journal of computational chemistry 2011, 32, 3423– 3432, DOI: 10.1002/jcc.21921
Google Scholar
48
Applying efficient implicit nongeometric constraints in alchemical free energy simulations
Knight, Jennifer L.; Brooks, Charles L.
Journal of Computational Chemistry (2011), 32 (16), 3423-3432CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in mol. dynamics simulations. Advanced variations of alchem. free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are assocd. with the progress variables in alchem. transformations are treated as dynamic variables and are constrained such that: 0 ≤ λi ≤ 1 and .sum.i λi = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by a functional form of λi that oscillates between λi = 0 and λi = 1 and has relatively steep transitions between these endpoints. Copyright for JCC Journal: © 2011 Wiley Periodicals, Inc.; J. Comput. Chem., 2011.
49
Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19– 25, DOI: 10.1016/j.softx.2015.06.001
Google Scholar
There is no corresponding record for this reference.
50
Singhal, A. K.; Chien, K. Y.; Wu, W. G.; Rule, G. S. Solution structure of cardiotoxin V from Naja naja atra. Biochemistry 1993, 32, 8036– 8044, DOI: 10.1021/bi00082a026
Google Scholar
50
Solution structure of cardiotoxin V from Naja naja atra
Singhal, Arun K.; Chien, Kun Yi; Wu, Wen Guey; Rule, Gordon S.
Biochemistry (1993), 32 (31), 8036-44CODEN: BICHAW; ISSN:0006-2960.
To det. whether the unique activity of this cardiotoxin (CTX) is attributable to its tertiary structure, the soln. structure of CTX V was detd. by NMR methods. On the basis of these studies, this cardiotoxin has the same general topol. as other members of the family, and thus its unusual properties do not arise from any gross structural differences that are detectable by soln. NMR methods. Mol. dynamics calcns. indicate that residues 36-50 show concerted fluctuations. On the basis of sequence similarity, the authors postulate that residues 30-34 are important in detg. the specificity of cardiotoxins for fusion vs. lysis of vesicles.
51
Ramanadham, M.; Sieker, L.; Jensen, L. Refinement of triclinic lysozyme: II. The method of stereochemically restrained least squares. Acta Crystallographica Section B: Structural Science 1990, 46, 63– 69, DOI: 10.1107/S0108768189009195
Google Scholar
51
Refinement of triclinic lysozyme: II. The method of stereochemically restrained least squares
Ramanadham M; Sieker L C; Jensen L H
Acta crystallographica. Section B, Structural science (1990), 46 ( Pt 1) (), 63-9 ISSN:0108-7681.
Refinement of triclinic lysozyme by restrained least squares against the 2 A resolution X-ray data is described, beginning with the model from cycle 17 of the preceding paper [Hodsdon, Brown, Sieker & Jensen (1990). Acta Cryst. B46, 54-62]. After 20 refinement cycles, R stood at 0.172. Nevertheless, serious errors involving both main-chain and side-chain atoms still remained, requiring numerous model rebuilding sessions interleaved with refinement cycles. After 63 cycles R = 0.124 for the model which includes all protein atoms, 249 water oxygen sites and five nitrate ions. Although the overall B is relatively low, 10.5 A2, B's for atoms in the region of residues 101-103, toward the termini of some of the longer side chains, and in the region of the C terminus of the main chain exceed 20 A2, indicating relatively high atomic mobilities, disorder, or remaining errors in the model.
52
Sauguet, L.; Poitevin, F.; Murail, S.; Van Renterghem, C.; Moraga-Cid, G.; Malherbe, L.; Thompson, A. W.; Koehl, P.; Corringer, P.-J.; Baaden, M.; Delarue, M. Structural basis for ion permeation mechanism in pentameric ligand-gated ion channels. EMBO journal 2013, 32, 728– 741, DOI: 10.1038/emboj.2013.17
Google Scholar
52
Structural basis for ion permeation mechanism in pentameric ligand-gated ion channels
Sauguet, Ludovic; Poitevin, Frederic; Murail, Samuel; Van Renterghem, Catherine; Moraga-Cid, Gustavo; Malherbe, Laurie; Thompson, Andrew W.; Koehl, Patrice; Corringer, Pierre-Jean; Baaden, Marc; Delarue, Marc
EMBO Journal (2013), 32 (5), 728-741CODEN: EMJODG; ISSN:0261-4189. (Nature Publishing Group)
To understand the mol. mechanism of ion permeation in pentameric ligand-gated ion channels (pLGIC), the structure of an open form of GLIC, a prokaryotic pLGIC, was solved at 2.4 Å. Anomalous diffraction data were used to place bound anions and cations. This reveals ordered water mols. at the level of two rings of hydroxylated residues (named Ser6' and Thr2') that contribute to the ion selectivity filter. Two water pentagons are obsd., a self-stabilized ice-like water pentagon and a second wider water pentagon, with one sodium ion between them. Single-channel electrophysiol. shows that the side-chain hydroxyl of Ser6' is crucial for ion translocation. Simulations and electrostatics calcns. complemented the description of hydration in the pore and suggest that the water pentagons obsd. in the crystal are important for the ion to cross hydrophobic constriction barriers. Simulations that pull a cation through the pore reveal that residue Ser6' actively contributes to ion translocation by reorienting its side chain when the ion is going through the pore. Generalization of these findings to the pLGIC family is proposed.
53
Lee, T.-W.; Qasim, M., Jr; Laskowski, M., Jr; James, M. N. Structural insights into the non-additivity effects in the sequence-to-reactivity algorithm for serine peptidases and their inhibitors. Journal of molecular biology 2007, 367, 527– 546, DOI: 10.1016/j.jmb.2007.01.008
Google Scholar
53
Structural Insights into the Non-additivity Effects in the Sequence-to-Reactivity Algorithm for Serine Peptidases and their Inhibitors
Lee, Ting-Wai; Qasim, M. A.; Laskowski, Michael; James, Michael N. G.
Journal of Molecular Biology (2007), 367 (2), 527-546CODEN: JMOBAK; ISSN:0022-2836. (Elsevier Ltd.)
Sequence-to-reactivity algorithms (SRAs) for proteins have the potential of being broadly applied in mol. design. Recently, Laskowski et al. have reported an additivity-based SRA that accurately predicts most of the std. free energy changes of assocn. for variants of turkey ovomucoid third domain (OMTKY3) with six serine peptidases, one of which is streptogrisin B (commonly known as Streptomyces griseus peptidase B, SGPB). Non-additivity effects for residues 18I and 32I, and for residues 20I and 32I of OMTKY3 occurred when the assocns. with SGPB were predicted using the SRA. To elucidate precisely the mechanics of these non-additivity effects in structural terms, we have detd. the crystal structures of the unbound OMTKY3 (with Gly32I as in the wild-type amino acid sequence) at a resoln. of 1.16 Å, the unbound Ala32I variant of OMTKY3 at a resoln. of 1.23 Å, and the SGPB:OMTKY3-Ala32I complex (equil. assocn. const. Ka = 7.1×109 M-1 at 21(±2) C°, pH 8.3) at a resoln. of 1.70 Å. Extensive comparisons with the crystal structure of the unbound OMTKY3 confirm our understanding of some previously addressed non-additivity effects. Unexpectedly, SGPB and OMTKY3-Ala32I form a 1:2 complex in the crystal. Comparison with the SGPB:OMTKY3 complex shows a conformational change in the SGPB:OMTKY3-Ala32I complex, resulting from a hinged rigid-body rotation of the inhibitor caused by the steric hindrance between the Me group of Ala32IA of the inhibitor and Pro192BE of the peptidase. This perturbs the interactions among residues 18I, 20I, 32I and 36I of the inhibitor, probably resulting in the above non-additivity effects. This conformational change also introduces residue 10I as an addnl. hyper-variable contact residue to the SRA.
54
Huang, J.; MacKerell, A. D., Jr CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data. Journal of computational chemistry 2013, 34, 2135– 2145, DOI: 10.1002/jcc.23354
Google Scholar
54
CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data
Huang, Jing; MacKerell, Alexander D. Jr
Journal of Computational Chemistry (2013), 34 (25), 2135-2145CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
Protein structure and dynamics can be characterized on the atomistic level with both NMR expts. and mol. dynamics (MD) simulations. Here, the authors quantify the ability of the recently presented CHARMM36 (C36) force field (FF) to reproduce various NMR observables using MD simulations. The studied NMR properties include backbone scalar couplings across hydrogen bonds, residual dipolar couplings (RDCs) and relaxation order parameter, as well as scalar couplings, RDCs, and order parameters for side-chain amino- and methyl-contg. groups. The C36 FF leads to better correlation with exptl. data compared to the CHARMM22/CMAP FF and suggest using C36 in protein simulations. Although both CHARMM FFs contains the same nonbond parameters, the authors' results show how the changes in the internal parameters assocd. with the peptide backbone via CMAP and the χ1 and χ2 dihedral parameters leads to improved treatment of the analyzed nonbond interactions. This highlights the importance of proper treatment of the internal covalent components in modeling nonbond interactions with mol. mechanics FFs.
55
Buslaev, P.; Aho, N.; Jansen, A.; Bauer, P.; Hess, B.; Groenhof, G. Best practices in constant pH MD simulations: accuracy and precision. J. Chem. Theory Comput. 2022, DOI: 10.1021/acs.jctc.2c00517 .
Google Scholar
There is no corresponding record for this reference.
56
Martini-website, General purpose coarse-grained force field. http://cgmartini.nl/index.php/tools2/proteins-and-bilayers/204-martinize (accessed October 01, 2021).
Google Scholar
There is no corresponding record for this reference.
57
Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J. Phys. Chem. B 2001, 105, 6474– 6487, DOI: 10.1021/jp003919d
Google Scholar
57
Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides
Kaminski, George A.; Friesner, Richard A.; Tirado-Rives, Julian; Jorgensen, William L.
Journal of Physical Chemistry B (2001), 105 (28), 6474-6487CODEN: JPCBFK; ISSN:1089-5647. (American Chemical Society)
We present results of improving the OPLS-AA force field for peptides by means of refitting the key Fourier torsional coeffs. The fitting technique combines using accurate ab initio data as the target, choosing an efficient fitting subspace of the whole potential-energy surface, and detg. wts. for each of the fitting points based on magnitudes of the potential-energy gradient. The av. energy RMS deviation from the LMP2/cc-pVTZ(-f)//HF/6-31G** data is reduced by ∼40% from 0.81 to 0.47 kcal/mol as a result of the fitting for the electrostatically uncharged dipeptides. Transferability of the parameters is demonstrated by using the same alanine dipeptide-fitted backbone torsional parameters for all of the other dipeptides (with the appropriate side-chain refitting) and the alanine tetrapeptide. Parameters of nonbonded interactions have also been refitted for the sulfur-contg. dipeptides (cysteine and methionine), and the validity of the new Coulombic charges and the van der Waals σ's and ε's is proved through reproducing gas-phase energies of complex formation heats of vaporization and densities of pure model liqs. Moreover, a novel approach to fitting torsional parameters for electrostatically charged mol. systems has been presented and successfully tested on five dipeptides with charged side chains.
58
Buck, M.; Bouguet-Bonnet, S.; Pastor, R. W.; MacKerell, A. D., Jr Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme. Biophysical journal 2006, 90, L36– L38, DOI: 10.1529/biophysj.105.078154
Google Scholar
58
Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme
Buck, Matthias; Bouguet-Bonnet, Sabine; Pastor, Richard W.; MacKerell, Alexander D., Jr.
Biophysical Journal (2006), 90 (4), L36-L38CODEN: BIOJAU; ISSN:0006-3495. (Biophysical Society)
The recently developed CMAP correction to the CHARMM22 force field (C22) is evaluated from 25 ns mol. dynamics simulations on hen lysozyme. Substantial deviations from exptl. backbone root mean-square fluctuations and N-H NMR order parameters obtained in the C22 trajectories (esp. in the loops) are eliminated by the CMAP correction. Thus, the C22/CMAP force field yields improved dynamical and structural properties of proteins in mol. dynamics simulations.
59
Durell, S. R.; Brooks, B. R.; Ben-Naim, A. Solvent-induced forces between two hydrophilic groups. J. Phys. Chem. 1994, 98, 2198– 2202, DOI: 10.1021/j100059a038
Google Scholar
59
Solvent-Induced Forces between Two Hydrophilic Groups
Durell, Stewart R.; Brooks, Bernard R.; Ben-Naim, Arieh
Journal of Physical Chemistry (1994), 98 (8), 2198-202CODEN: JPCHAX; ISSN:0022-3654.
Mol. dynamics simulations were used to calc. the force between two simple hydrophilic solutes in dil. aq. soln. The "solutes" were two water mols. in the same relative orientation as the next-nearest neighbors in hexagonal ice I. Both the direct and solvent-induced contributions to the force were calcd. as a function of sepn. distance. The total force between the solutes was found to be most attractive at 5.0 Å (-1.6 kcal/mol/Å). The potential of mean force had a min. at 4.3 Å, which is 0.2 Å closer than the next-nearest-neighbor distance in ice. A parallel set of simulations were conducted with the partial charges on the "solutes" removed to examine hydrophobic analogs. In this case, the total force was most attractive at 3.5 Å (-0.9 kcal/mol/Å), and the min. of the potential was at the contact distance of 3.2 Å. In agreement with earlier predictions, the max. solvent-induced contribution to the potential was ca. 4 times more neg. for the hydrophilic "solutes" than for the hydrophobic ones. These differences are shown to be due predominantly to a solvent water mol. which simultaneously hydrogen bonds to both hydrophilic "solutes". The results support earlier assertions that solvent-induced interactions between polar amino acid residues are more important in protein folding and stability than generally considered.
60
Neria, E.; Fischer, S.; Karplus, M. Simulation of activation free energies in molecular systems. J. Chem. Phys. 1996, 105, 1902– 1921, DOI: 10.1063/1.472061
Google Scholar
60
Simulation of activation free energies in molecular systems
Neria, Eyal; Fischer, Stefan; Karplus, Martin
Journal of Chemical Physics (1996), 105 (5), 1902-1921CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
A method is presented for detg. activation free energies in complex mol. systems. The method relies on knowledge of the min. energy path and bases the activation free energy calcn. on moving along this path from a min. to a saddle point. Use is made of a local reaction coordinate which describes the advance of the reaction in each segment of the min. energy path. The activation free energy is formulated as a sum of two terms. The first is due to the change in the local reaction coordinate between the endpoints of each segment of the path. The second is due to the change in direction of the min. energy path between consecutive segments. Both contributions can be obtained by mol. dynamics simulations with a constraint on the local reaction coordinate. The method is illustrated by applying it to a model potential and to the C7eq to C7ax transition in the alanine dipeptide. It is found that the term due to the change of direction in the reaction path can make a substantial contribution to the activation free energy.
61
Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable water model for the coarse-grained MARTINI force field. PLoS computational biology 2010, 6, e1000810, DOI: 10.1371/journal.pcbi.1000810
Google Scholar
61
Polarizable water model for the coarse-grained MARTINI force field
Yesylevskyy Semen O; Schafer Lars V; Sengupta Durba; Marrink Siewert J
PLoS computational biology (2010), 6 (6), e1000810 ISSN:.
Coarse-grained (CG) simulations have become an essential tool to study a large variety of biomolecular processes, exploring temporal and spatial scales inaccessible to traditional models of atomistic resolution. One of the major simplifications of CG models is the representation of the solvent, which is either implicit or modeled explicitly as a van der Waals particle. The effect of polarization, and thus a proper screening of interactions depending on the local environment, is absent. Given the important role of water as a ubiquitous solvent in biological systems, its treatment is crucial to the properties derived from simulation studies. Here, we parameterize a polarizable coarse-grained water model to be used in combination with the CG MARTINI force field. Using a three-bead model to represent four water molecules, we show that the orientational polarizability of real water can be effectively accounted for. This has the consequence that the dielectric screening of bulk water is reproduced. At the same time, we parameterized our new water model such that bulk water density and oil/water partitioning data remain at the same level of accuracy as for the standard MARTINI force field. We apply the new model to two cases for which current CG force fields are inadequate. First, we address the transport of ions across a lipid membrane. The computed potential of mean force shows that the ions now naturally feel the change in dielectric medium when moving from the high dielectric aqueous phase toward the low dielectric membrane interior. In the second application we consider the electroporation process of both an oil slab and a lipid bilayer. The electrostatic field drives the formation of water filled pores in both cases, following a similar mechanism as seen with atomistically detailed models.
62
Berendsen, H. J.; Postma, J. P.; van Gunsteren, W. F.; Hermans, J. Intermolecular forces; Springer, 1981; pp 331– 342.
Google Scholar
There is no corresponding record for this reference.
63
De Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini straight: Boosting performance using a shorter cutoff and GPUs. Comput. Phys. Commun. 2016, 199, 1– 7, DOI: 10.1016/j.cpc.2015.09.014
Google Scholar
63
Martini straight: Boosting performance using a shorter cutoff and GPUs
de Jong, Djurre H.; Baoukina, Svetlana; Ingolfsson, Helgi I.; Marrink, Siewert J.
Computer Physics Communications (2016), 199 (), 1-7CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
In mol. dynamics simulations, sufficient sampling is of key importance and a continuous challenge in the field. The coarse grain Martini force field has been widely used to enhance sampling. In its original implementation, this force field applied a shifted Lennard-Jones potential for the non-bonded van der Waals interactions, to avoid problems related to a relatively short cutoff. Here we investigate the use of a straight cutoff Lennard-Jones potential with potential modifiers. Together with a Verlet neighbor search algorithm, the modified potential allows the use of GPUs to accelerate the computations in Gromacs. We find that this alternative potential has little influence on most of the properties studied, including partitioning free energies, bulk liq. properties and bilayer properties. At the same time, energy conservation is kept within reasonable bounds. We conclude that the newly proposed straight cutoff approach is a viable alternative to the std. shifted potentials used in Martini, offering significant speedup even in the absence of GPUs.
64
Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101, DOI: 10.1063/1.2408420
Google Scholar
64
Canonical sampling through velocity rescaling
Bussi, Giovanni; Donadio, Davide; Parrinello, Michele
Journal of Chemical Physics (2007), 126 (1), 014101/1-014101/7CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
The authors present a new mol. dynamics algorithm for sampling the canonical distribution. In this approach the velocities of all the particles are rescaled by a properly chosen random factor. The algorithm is formally justified and it is shown that, in spite of its stochastic nature, a quantity can still be defined that remains const. during the evolution. In numerical applications this quantity can be used to measure the accuracy of the sampling. The authors illustrate the properties of this new method on Lennard-Jones and TIP4P water models in the solid and liq. phases. Its performance is excellent and largely independent of the thermostat parameter also with regard to the dynamic properties.
65
Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182– 7190, DOI: 10.1063/1.328693
Google Scholar
65
Polymorphic transitions in single crystals: A new molecular dynamics method
Parrinello, M.; Rahman, A.
Journal of Applied Physics (1981), 52 (12), 7182-90CODEN: JAPIAU; ISSN:0021-8979.
A Lagrangian formulation is introduced; it can be used to make mol. dynamics (MD) calcns. on systems under the most general, externally applied, conditions of stress. In this formulation the MD cell shape and size can change according to dynamic equations given by this Lagrangian. This MD technique was used to the study of structural transitions of a Ni single crystal under uniform uniaxial compressive and tensile loads. Some results regarding the stress-strain relation obtained by static calcns. are invalid at finite temp. Under compressive loading, the model of Ni shows a bifurcation in its stress-strain relation; this bifurcation provides a link in configuration space between cubic and hexagonal close packing. Such a transition could perhaps be obsd. exptl. under extreme conditions of shock.
66
Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: a linear constraint solver for molecular simulations. Journal of computational chemistry 1997, 18, 1463– 1472, DOI: 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H
Google Scholar
66
LINCS: a linear constraint solver for molecular simulations
Hess, Berk; Bekker, Henk; Berendsen, Herman J. C.; Fraaije, Johannes G. E. M.
Journal of Computational Chemistry (1997), 18 (12), 1463-1472CODEN: JCCHDD; ISSN:0192-8651. (Wiley)
We present a new LINear Constraint Solver (LINCS) for mol. simulations with bond constraints using the enzyme lysozyme and a 32-residue peptide as test systems. The algorithm is inherently stable, as the constraints themselves are reset instead of derivs. of the constraints, thereby eliminating drift. Although the derivation of the algorithm is presented in terms of matrixes, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large mols. At the same accuracy, the LINCS algorithm is 3-4 times faster than the SHAKE algorithm. Parallelization of the algorithm is straightforward.
67
Miyamoto, S.; Kollman, P. A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. Journal of computational chemistry 1992, 13, 952– 962, DOI: 10.1002/jcc.540130805
Google Scholar
67
SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models
Miyamoto, Shuichi; Kollman, Peter A.
Journal of Computational Chemistry (1992), 13 (8), 952-62CODEN: JCCHDD; ISSN:0192-8651.
An anal. algorithm, called SETTLE, for resetting the positions and velocities to satisfy the holonomic constraints on the rigid water model is presented. This method is based on the Cartesian coordinate system and can be used in place of SHAKE and RATTLE. The authors implemented this algorithm in the SPASMS package of mol. mechanics and dynamics. Several series of mol. dynamics simulations were carried out to examine the performance of the new algorithm in comparison with the original RATTLE method. SETTLE is of higher accuracy and is faster than RATTLE with reasonable tolerances by three to nine times on a scalar machine. The performance improvement ranged from factors of 26 to 98 on a vector machine since the method presented is not iterative.
68
Hub, J. S.; de Groot, B. L.; Grubmüller, H.; Groenhof, G. Quantifying artifacts in Ewald simulations of inhomogeneous systems with a net charge. J. Chem. Theory Comput. 2014, 10, 381– 390, DOI: 10.1021/ct400626b
Google Scholar
68
Quantifying Artifacts in Ewald Simulations of Inhomogeneous Systems with a Net Charge
Hub, Jochen S.; de Groot, Bert L.; Grubmueller, Helmut; Groenhof, Gerrit
Journal of Chemical Theory and Computation (2014), 10 (1), 381-390CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
Ewald summation, which has become the de facto std. for computing electrostatic interactions in biomol. simulations, formally requires that the simulation box is neutral. For non-neutral systems, the Ewald algorithm implicitly introduces a uniform background charge distribution that effectively neutralizes the simulation box. Because a uniform distribution of counter charges typically deviates from the spatial distribution of counterions in real systems, artifacts may arise, in particular in systems with an inhomogeneous dielec. const. Here, we derive an anal. expression for the effect of using an implicit background charge instead of explicit counterions, on the chem. potential of ions in heterogeneous systems, which (i) provides a quant. criterion for deciding if the background charge offers an acceptable trade-off between artifacts arising from sampling problems and artifacts arising from the homogeneous background charge distribution, and (ii) can be used to correct this artifact in certain cases. Our model quantifies the artifact in terms of the difference in charge d. between the non-neutral system with a uniform neutralizing background charge and the real neutral system with a phys. correct distribution of explicit counterions. We show that for inhomogeneous systems, such as proteins and membranes in water, the artifact manifests itself by an overstabilization of ions inside the lower dielec. by tens to even hundreds kilojoules per mol. We have tested the accuracy of our model in mol. dynamics simulations and found that the error in the calcd. free energy for moving a test charge from water into hexadecane, at different net charges of the system and different simulation box sizes, is correctly predicted by the model. The calcns. further confirm that the incorrect distribution of counter charges in the simulation box is solely responsible for the errors in the PMFs.
69
Thurlkill, R. L.; Grimsley, G. R.; Scholtz, J. M.; Pace, C. N. pK values of the ionizable groups of proteins. Protein science 2006, 15, 1214– 1218, DOI: 10.1110/ps.051840806
Google Scholar
69
pK values of the ionizable groups of proteins
Thurlkill, Richard L.; Grimsley, Gerald R.; Scholtz, J. Martin; Pace, C. Nick
Protein Science (2006), 15 (5), 1214-1218CODEN: PRCIEI; ISSN:0961-8368. (Cold Spring Harbor Laboratory Press)
We have used potentiometric titrns. to measure the pK values of the ionizable groups of proteins in alanine pentapeptides with appropriately blocked termini. These pentapeptides provide an improved model for the pK values of the ionizable groups in proteins. Our pK values detd. in 0.1 M KCl at 25°C are: 3.67 ± 0.03 (α-carboxyl), 3.67 ± 0.04 (Asp), 4.25 ± 0.05 (Glu), 6.54 ± 0.04 (His), 8.00 ± 0.03 (α-amino), 8.55 ± 0.03 (Cys), 9.84 ± 0.11 (Tyr), and 10.40 ± 0.08 (Lys). The pK values of some groups differ from the Nozaki and Tanford (N&T) pK values often used in the literature: Asp (3.67 this work vs. 4.0 N&T); His (6.54 this work vs. 6.3 N&T); α-amino (8.00 this work vs. 7.5 N&T); Cys (8.55 this work vs. 9.5 N&T); and Tyr (9.84 this work vs. 9.6 N&T). Our pK values will be useful to those who study pK perturbations in folded and unfolded proteins, and to those who use theory to gain a better understanding of the factors that det. the pK values of the ionizable groups of proteins.
70
Tanokura, M. 1H-NMR study on the tautomerism of the imidazole ring of histidine residues: I. Microscopic pK values and molar ratios of tautomers in histidine-containing peptides. Biochimica et Biophysica Acta (BBA)-Protein Structure and Molecular Enzymology 1983, 742, 576– 585, DOI: 10.1016/0167-4838(83)90276-5
Google Scholar
There is no corresponding record for this reference.
71
Chiang, C.-M.; Chien, K.-Y.; Lin, H.-j.; Lin, J.-F.; Yeh, H.-C.; Ho, P.-l.; Wu, W.-g. Conformational change and inactivation of membrane phospholipid-related activity of cardiotoxin V from Taiwan cobra venom at acidic pH. Biochemistry 1996, 35, 9167– 9176, DOI: 10.1021/bi952823k
Google Scholar
71
Conformational Change and Inactivation of Membrane Phospholipid-Related Activity of Cardiotoxin V from Taiwan Cobra Venom at Acidic pH
Chiang, Chien-Min; Chien, Kun-Yi; Lin, Hai-jui; Lin, Ji-Fu; Yeh, Hsien-Chi; Ho, Pei-li; Wu, Wen-guey
Biochemistry (1996), 35 (28), 9167-9176CODEN: BICHAW; ISSN:0006-2960. (American Chemical Society)
The phospholipid binding activity of cardiotoxin V from Naja naja atra (CTX A5) was studied by use of Langmuir monolayers and found to exhibit pH-dependence in binding to phosphatidylcholine membrane with an apparent pKa around 6.0. Proton NMR investigation of the CTX A5 mol. in the presence of phosphatidylcholine micelles reveals a decrease in assocn. of CTX A5 with membranes at low pH as a result of the protonation of His-4 near the membrane binding site of loop I region of CTX. The pH-dependent binding can be attributed mainly, but not solely, to the change in charge content of the CTX mol. upon His-4 protonation at the membrane/water interface. This is shown by analyzing the pH- and ionic strength dependence of binding of CTXs to phospholipid monolayers according to Gouy-Chapman theory. The protonation of the His-4 residue also results in a local conformational change in the loop I region since the chem. shifts of amide protons for the amino acid residues from Cys-3 to Thr-14 are all found to vary as a function of pH with an apparent pKa similar to that of His-4. Interestingly, the effect is relayed to other amino acid residues in the structural core of the protein such as those in C-terminal (Lys-60, Cys-61, and Asn-62) and triple-stranded antiparallel β-sheet (Cys-22, Lys-24, Ala-25, Arg-38, and Ala-41) regions. An addnl. local conformational change in the mol. results around pH 5 as evidenced by CD spectroscopic studies, although this change does not affect the characteristic β-sheet and three-finger loop structure of CTX mol. as revealed by two-dimensional NOESY 1H NMR study. The latter conformational change at acidic pH, however, completely inactivates CTX-induced aggregation/fusion activity of sphingomyelin vesicles. The results suggest that deciphering the functional sites of CTXs on the basis of structure and dynamics detd. at low pH should be done with caution. Since 19 out of 44 CTX homologs with known amino acid sequence contain His-4, the effect of His-4 on the structure and function of CTX mols. is important and is discussed in terms of the diverse membrane targets of CTX subtypes. Also discussed is the pH-induced activation of snake venom proteins in the victim.
72
Chiang, C.-M.; Chang, S.-L.; Lin, H.-j.; Wu, W.-g. The role of acidic amino acid residues in the structural stability of snake cardiotoxins. Biochemistry 1996, 35, 9177– 9186, DOI: 10.1021/bi960077t
Google Scholar
72
The Role of Acidic Amino Acid Residues in the Structural Stability of Snake Cardiotoxins
Chiang, Chien-Min; Chang, Shou-Lin; Lin, Hai-jui; Wu, Wen-guey
Biochemistry (1996), 35 (28), 9177-9186CODEN: BICHAW; ISSN:0006-2960. (American Chemical Society)
PH-induced structural transitions occurred in cardiotoxin V from Naja naja atra (CTX A5) at two pH values as judged by the CD ellipticity around 195 nm: an increase in the β-sheet content occurred around pH 4 and followed by a decrease, therein, around pH 2. The pKa of three acidic amino acid residues in CTX A5, i.e., Glu-17, Asp-42, and Asp-59, were detd. to be 4.0, 3.2, and below 2.3, resp., by NMR spectroscopy. The low pKa value of Asp-59 implies salt bridge formation between Lys-2 and Asp-59. Thus, electrostatic interaction may stabilize the three loop structure in addn. to the hydrogen bonds between N- and C-termini of CTX mol. Second, 2,2,2-trifluoroethanol (TFE) and guanidinium chloride (GdmHCl) were found to induce α-helical and random coil formation, resp., in CTX A5 and eight other β-sheet CTXs. Comparison of the relative potencies of TFE and GdmHCl to induce structural changes suggests that the amino acid residue located at position 17 plays a role in the structural stability. Specifically, CTXs contg. neg. charged Glu-17 are least stable. It is suggested that Glu-17 may perturb the interaction between Lys-2 and Asp-59, and thus the overall stability of β-sheet, in the presence of denaturing reagent. In conclusion, the perturbed structural stability of CTXs may partially explain the lower activity CTX exhibits at acidic pH. A structural model to account for the unfolding and refolding of CTX mols. without the breaking of disulfide bonds is also proposed.
73
Webb, H.; Tynan-Connolly, B. M.; Lee, G. M.; Farrell, D.; O’Meara, F.; Søndergaard, C. R.; Teilum, K.; Hewage, C.; McIntosh, L. P.; Nielsen, J. E. Remeasuring HEWL pKa values by NMR spectroscopy: Methods, analysis, accuracy, and implications for theoretical pKa calculations. Proteins: Struct., Funct., Bioinf. 2011, 79, 685– 702, DOI: 10.1002/prot.22886
Google Scholar
73
Remeasuring HEWL pKa values by NMR spectroscopy: methods, analysis, accuracy, and implications for theoretical pKa calculations
Webb, Helen; Tynan-Connolly, Barbara Mary; Lee, Gregory M.; Farrell, Damien; O'Meara, Fergal; Sondergaard, Chresten R.; Teilum, Kaare; Hewage, Chandralal; McIntosh, Lawrence P.; Nielsen, Jens Erik
Proteins: Structure, Function, and Bioinformatics (2011), 79 (3), 685-702CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
Site-specific pKa values measured by NMR spectroscopy provide essential information on protein electrostatics, the pH-dependence of protein structure, dynamics and function, and constitute an important benchmark for protein pKa calcn. algorithms. Titrn. curves can be measured by tracking the NMR chem. shifts of several reporter nuclei vs. sample pH. However, careful anal. of these curves is needed to ext. residue-specific pKa values since pH-dependent chem. shift changes can arise from many sources, including through-bond inductive effects, through-space elec. field effects, and conformational changes. We have re-measured titrn. curves for all carboxylates and His 15 in Hen Egg White Lysozyme (HEWL) by recording the pH-dependent chem. shifts of all backbone amide nitrogens and protons, Asp/Glu side chain protons and carboxyl carbons, and imidazole protonated carbons and protons in this protein. We extd. pKa values from the resulting titrn. curves using std. fitting methods, and compared these values to each other, and with those measured previously by 1H NMR. This anal. gives insights into the true accuracy assocd. with exptl. measured pKa values. We find that apparent pKa values frequently differ by 0.5-1.0 units depending upon the nuclei monitored, and that larger differences occasionally can be obsd. The variation in measured pKa values, which reflects the difficulty in fitting and assigning pH-dependent chem. shifts to specific ionization equil., has significant implications for the exptl. procedures used for measuring protein pKa values, for the benchmarking of protein pKa calcn. algorithms, and for the understanding of protein electro-statics in general.
74
Lee, J.; Miller, B. T.; Damjanovic, A.; Brooks, B. R. Enhancing constant-pH simulation in explicit solvent with a two-dimensional replica exchange method. J. Chem. Theory Comput. 2015, 11, 2560– 2574, DOI: 10.1021/ct501101f
Google Scholar
74
Enhancing Constant-pH Simulation in Explicit Solvent with a Two-Dimensional Replica Exchange Method
Lee, Juyong; Miller, Benjamin T.; Damjanovic, Ana; Brooks, Bernard R.
Journal of Chemical Theory and Computation (2015), 11 (6), 2560-2574CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
We present a new method for enhanced sampling for const.-pH simulations in explicit water based on a two-dimensional (2D) replica exchange scheme. The new method is a significant extension of our previously developed const.-pH simulation method, which is based on enveloping distribution sampling (EDS) coupled with a one-dimensional (1D) Hamiltonian exchange method (HREM). EDS constructs a hybrid Hamiltonian from multiple discrete end state Hamiltonians that, in this case, represent different protonation states of the system. The ruggedness and heights of the hybrid Hamiltonian's energy barriers can be tuned by the smoothness parameter. Within the context of the 1D EDS-HREM method, exchanges are performed between replicas with different smoothness parameters, allowing frequent protonation-state transitions and sampling of conformations that are favored by the end-state Hamiltonians. In this work, the 1D method is extended to 2D with an addnl. dimension, external pH. Within the context of the 2D method (2D EDS-HREM), exchanges are performed on a lattice of Hamiltonians with different pH conditions and smoothness parameters. We demonstrate that both the 1D and 2D methods exactly reproduce the thermodn. properties of the semigrand canonical (SGC) ensemble of a system at a given pH. We have tested our new 2D method on aspartic acid, glutamic acid, lysine, a four residue peptide (sequence KAAE), and snake cardiotoxin. In all cases, the 2D method converges faster and without loss of precision; the only limitation is a loss of flexibility in how CPU time is employed. The results for snake cardiotoxin demonstrate that the 2D method enhances protonation-state transitions, samples a wider conformational space with the same amt. of computational resources, and converges significantly faster overall than the original 1D method.
75
Pezeshkian, W.; König, M.; Wassenaar, T. A.; Marrink, S. J. Backmapping triangulated surfaces to coarse-grained membrane models. Nat. Commun. 2020, 11, 2296, DOI: 10.1038/s41467-020-16094-y
Google Scholar
75
Backmapping triangulated surfaces to coarse-grained membrane models
Pezeshkian, Weria; Koenig, Melanie; Wassenaar, Tsjerk A.; Marrink, Siewert J.
Nature Communications (2020), 11 (1), 2296CODEN: NCAOBW; ISSN:2041-1723. (Nature Research)
Many biol. processes involve large-scale changes in membrane shape. Computer simulations of these processes are challenging since they occur across a wide range of spatiotemporal scales that cannot be investigated in full by any single current simulation technique. A potential soln. is to combine different levels of resoln. through a multiscale scheme. Here, we present a multiscale algorithm that backmaps a continuum membrane model represented as a dynamically triangulated surface (DTS) to its corresponding mol. model based on the coarse-grained (CG) Martini force field. Thus, we can use DTS simulations to equilibrate slow large-scale membrane conformational changes and then explore the local properties at CG resoln. We demonstrate the power of our method by backmapping a vesicular bud induced by binding of Shiga toxin and by transforming the membranes of an entire mitochondrion to near-at. resoln. Our approach opens the way to whole cell simulations at mol. detail.
76
Wolfram Research, Inc. Mathematica, Version 11.3. Champaign: IL, 2018.
Google Scholar
There is no corresponding record for this reference.

Cited By

Click to copy section linkSection link copied!

Explore this article's citation statements on scite.ai

This article is cited by 75 publications.

Tomás F. D. Silva, Giovanni Bussi. Characterizing RNA Oligomers Using Stochastic Titration Constant-pH Metadynamics Simulations. Journal of Chemical Information and Modeling 2025, 65 (7) , 3568-3580. https://doi.org/10.1021/acs.jcim.4c02185
Shijie Xu, Akira Onoda. Accurate and Rapid Prediction of Protein pKa: Protein Language Models Reveal the Sequence–pKa Relationship. Journal of Chemical Theory and Computation 2025, 21 (7) , 3752-3764. https://doi.org/10.1021/acs.jctc.4c01288
Mohannad J. Yousef, Nuno F. B. Oliveira, João N. M. Vitorino, Pedro B. P. S. Reis, Piotr Draczkowski, Maciej Maj, Krzysztof Jozwiak, Miguel Machuqueiro. Toward Accurate pH-Dependent Binding Constant Predictions Using Molecular Docking and Constant-pH MD Calculations. Journal of Chemical Theory and Computation 2025, 21 (5) , 2655-2667. https://doi.org/10.1021/acs.jctc.4c01291
Markéta Paloncýová, Mariana Valério, Ricardo Nascimento Dos Santos, Petra Kührová, Martin Šrejber, Petra Čechová, Dimitar A. Dobchev, Akshay Balsubramani, Pavel Banáš, Vikram Agarwal, Paulo C. T. Souza, Michal Otyepka. Computational Methods for Modeling Lipid-Mediated Active Pharmaceutical Ingredient Delivery. Molecular Pharmaceutics 2025, 22 (3) , 1110-1141. https://doi.org/10.1021/acs.molpharmaceut.4c00744
Eliane Briand, Bartosz Kohnke, Carsten Kutzner, Helmut Grubmüller. Constant pH Simulation with FMM Electrostatics in GROMACS. (A) Design and Applications. Journal of Chemical Theory and Computation 2025, 21 (4) , 1762-1786. https://doi.org/10.1021/acs.jctc.4c01318
Bartosz Kohnke, Eliane Briand, Carsten Kutzner, Helmut Grubmüller. Constant pH Simulation with FMM Electrostatics in GROMACS. (B) GPU Accelerated Hamiltonian Interpolation. Journal of Chemical Theory and Computation 2025, 21 (4) , 1787-1804. https://doi.org/10.1021/acs.jctc.4c01319
David Beyer, Pablo M. Blanco, Jonas Landsgesell, Peter Košovan, Christian Holm. How To Correct Erroneous Symmetry-Breaking in Coarse-Grained Constant-pH Simulations. Journal of Chemical Theory and Computation 2025, 21 (3) , 1396-1404. https://doi.org/10.1021/acs.jctc.4c01010
Marc Domingo, Horacio V. Guzman, Matej Kanduč, Jordi Faraudo. Electrostatic Interaction between SARS-CoV-2 and Charged Surfaces: Spike Protein Evolution Changed the Game. Journal of Chemical Information and Modeling 2025, 65 (1) , 240-251. https://doi.org/10.1021/acs.jcim.4c01724
Craig A. Peeples, Ruibin Liu, Jana Shen. Force Field Limitations of All-Atom Continuous Constant pH Molecular Dynamics. The Journal of Physical Chemistry B 2024, 128 (47) , 11616-11624. https://doi.org/10.1021/acs.jpcb.4c05971
Richard S. Hong, Busayo D. Alagbe, Alessandra Mattei, Ahmad Y. Sheikh, Mark E. Tuckerman. Enhanced and Efficient Predictions of Dynamic Ionization through Constant-pH Adiabatic Free Energy Dynamics. Journal of Chemical Theory and Computation 2024, 20 (22) , 10010-10021. https://doi.org/10.1021/acs.jctc.4c00704
Noora Aho, Gerrit Groenhof, Pavel Buslaev. What Is the Protonation State of Proteins in Crystals? Insights from Constant pH Molecular Dynamics Simulations. The Journal of Physical Chemistry B 2024, 128 (45) , 11124-11133. https://doi.org/10.1021/acs.jpcb.4c05947
Zhongxin Cui, Yuefeng Wang, Lei Zhang, Haishan Qi. Zwitterionic Peptides: From Mechanism, Design Strategies to Applications. ACS Applied Materials & Interfaces 2024, 16 (42) , 56497-56518. https://doi.org/10.1021/acsami.4c08891
Wonmuk Hwang, Steven L. Austin, Arnaud Blondel, Eric D. Boittier, Stefan Boresch, Matthias Buck, Joshua Buckner, Amedeo Caflisch, Hao-Ting Chang, Xi Cheng, Yeol Kyo Choi, Jhih-Wei Chu, Michael F. Crowley, Qiang Cui, Ana Damjanovic, Yuqing Deng, Mike Devereux, Xinqiang Ding, Michael F. Feig, Jiali Gao, David R. Glowacki, James E. Gonzales, II, Mehdi Bagerhi Hamaneh, Edward D. Harder, Ryan L. Hayes, Jing Huang, Yandong Huang, Phillip S. Hudson, Wonpil Im, Shahidul M. Islam, Wei Jiang, Michael R. Jones, Silvan Käser, Fiona L. Kearns, Nathan R. Kern, Jeffery B. Klauda, Themis Lazaridis, Jinhyuk Lee, Justin A. Lemkul, Xiaorong Liu, Yun Luo, Alexander D. MacKerell, Jr., Dan T. Major, Markus Meuwly, Kwangho Nam, Lennart Nilsson, Victor Ovchinnikov, Emanuele Paci, Soohyung Park, Richard W. Pastor, Amanda R. Pittman, Carol Beth Post, Samarjeet Prasad, Jingzhi Pu, Yifei Qi, Thenmalarchelvi Rathinavelan, Daniel R. Roe, Benoit Roux, Christopher N. Rowley, Jana Shen, Andrew C. Simmonett, Alexander J. Sodt, Kai Töpfer, Meenu Upadhyay, Arjan van der Vaart, Luis Itza Vazquez-Salazar, Richard M. Venable, Luke C. Warrensford, H. Lee Woodcock, Yujin Wu, Charles L. Brooks, III, Bernard R. Brooks, Martin Karplus. CHARMM at 45: Enhancements in Accessibility, Functionality, and Speed. The Journal of Physical Chemistry B 2024, 128 (41) , 9976-10042. https://doi.org/10.1021/acs.jpcb.4c04100
Junjie Shen, Yi Sun, Xuanzhe Liu, Yimin Chai, Chunyang Wang, Jia Xu. Nerve Regeneration Potential of Antioxidant-Modified Black Phosphorus Quantum Dots in Peripheral Nerve Injury. ACS Nano 2024, 18 (34) , 23518-23536. https://doi.org/10.1021/acsnano.4c07285
Tomasz K. Piskorz, Laura Perez-Chirinos, Baofu Qiao, Ivan R. Sasselli. Tips and Tricks in the Modeling of Supramolecular Peptide Assemblies. ACS Omega 2024, 9 (29) , 31254-31273. https://doi.org/10.1021/acsomega.4c02628
Magnus Lundborg, Christian Wennberg, Erik Lindahl, Lars Norlén. Simulating the Skin Permeation Process of Ionizable Molecules. Journal of Chemical Information and Modeling 2024, 64 (13) , 5295-5302. https://doi.org/10.1021/acs.jcim.4c00722
Jonathan Lasham, Amina Djurabekova, Volker Zickermann, Janet Vonck, Vivek Sharma. Role of Protonation States in the Stability of Molecular Dynamics Simulations of High-Resolution Membrane Protein Structures. The Journal of Physical Chemistry B 2024, 128 (10) , 2304-2316. https://doi.org/10.1021/acs.jpcb.3c07421
Anton Jansen, Noora Aho, Gerrit Groenhof, Pavel Buslaev, Berk Hess. phbuilder: A Tool for Efficiently Setting up Constant pH Molecular Dynamics Simulations in GROMACS. Journal of Chemical Information and Modeling 2024, 64 (3) , 567-574. https://doi.org/10.1021/acs.jcim.3c01313
Tao Luo, Sébastien Le Crom, N. Tan Luong, Khalil Hanna, Jean-François Boily. Goethite-Bound Copper Controls the Fate of Antibiotics in Aquatic Environments. ACS ES&T Water 2024, 4 (2) , 638-647. https://doi.org/10.1021/acsestwater.3c00666
Luise Jacobsen, Laura Lydersen, Himanshu Khandelia. ATP-Bound State of the Uncoupling Protein 1 (UCP1) from Molecular Simulations. The Journal of Physical Chemistry B 2023, 127 (45) , 9685-9696. https://doi.org/10.1021/acs.jpcb.3c03473
Luís Borges-Araújo, Ilias Patmanidis, Akhil P. Singh, Lucianna H. S. Santos, Adam K. Sieradzan, Stefano Vanni, Cezary Czaplewski, Sergio Pantano, Wataru Shinoda, Luca Monticelli, Adam Liwo, Siewert J. Marrink, Paulo C. T. Souza. Pragmatic Coarse-Graining of Proteins: Models and Applications. Journal of Chemical Theory and Computation 2023, 19 (20) , 7112-7135. https://doi.org/10.1021/acs.jctc.3c00733
Varun Mandalaparthy, Madhusmita Tripathy, Nico F. A. van der Vegt. Anions and Cations Affect Amino Acid Dissociation Equilibria via Distinct Mechanisms. The Journal of Physical Chemistry Letters 2023, 14 (41) , 9250-9256. https://doi.org/10.1021/acs.jpclett.3c02062
Mykyta V. Prud, Alexander Kyrychenko, Oleg N. Kalugin. pH-Controllable Coating of Silver Nanoparticles with PMMA-b-PDMAEMA Oligomers: A Molecular Dynamics Simulation Study. The Journal of Physical Chemistry C 2023, 127 (24) , 11748-11759. https://doi.org/10.1021/acs.jpcc.3c02779
Zhitao Cai, Tengzi Liu, Qiaoling Lin, Jiahao He, Xiaowei Lei, Fangfang Luo, Yandong Huang. Basis for Accurate Protein pKa Prediction with Machine Learning. Journal of Chemical Information and Modeling 2023, 63 (10) , 2936-2947. https://doi.org/10.1021/acs.jcim.3c00254
Pavel Buslaev, Noora Aho, Anton Jansen, Paul Bauer, Berk Hess, Gerrit Groenhof. Best Practices in Constant pH MD Simulations: Accuracy and Sampling. Journal of Chemical Theory and Computation 2022, 18 (10) , 6134-6147. https://doi.org/10.1021/acs.jctc.2c00517
Jacilene Silva, Matheus Nunes da Rocha, Victor Moreira de Oliveira, Caio Henrique Alexandre Roberto, Francisco Nithael Melo Lúcio, Márcia Machado Marinho, Hélcio Silva dos Santos, Emmanuel Silva Marinho. Allosteric modulation of laeviganoid-based clerodane diterpenes derivatives in muscarinic acetylcholine M1 receptor against tinnitus: a structure-based virtual screening approach. Future Journal of Pharmaceutical Sciences 2025, 11 (1) https://doi.org/10.1186/s43094-025-00783-w
Yingqi Qiu, Jiahao Lu, Chenhao Zhao, Yuqiang Xiang, Aiqun Wu, Liqun Shen, Haiou Jiang. Synthesis of N-acylsulfonamide chromone derivatives as efficient anti-Candida albicans agents. Journal of Molecular Structure 2025, 1334 , 141887. https://doi.org/10.1016/j.molstruc.2025.141887
Aydin J. Hodala, Ian G. Wood, Paola Carbone. Understanding the influence of electrostatic interactions on observed pK shifts in surfactant aggregates using classical simulations. Journal of Molecular Liquids 2025, 427 , 127410. https://doi.org/10.1016/j.molliq.2025.127410
C. Risueño, I. Carbajo, D. Charro, N. G. A. Abrescia, I. Coluzza. pH-antenna residues trigger a large-scale conformational change in the large extracellular loop domain of the CD81 human receptor. The Journal of Chemical Physics 2025, 162 (15) https://doi.org/10.1063/5.0251361
Haiyan Xiong, Qiyuan Chang, Jiayi Ding, Shuyuan Wang, Wenhao Zhang, Yu Li, Yaochen Wu, Pengyan Lin, Chengyu Yang, Miaoxing Liu, Guicun Fang, Yiwei Yang, Jiongfang Xie, Dong Qi, Tao Jiang, Wenfeng Fu, Fen Hu, Yiming Chen, Rongcai Yue, Yanbin Li, Yong Cui, Min Li, Shilong Fan, Yufeng Yang, Yunlu Xu, Dong Li, Fenghua Zhang, Hu Zhao, Congxian Wu, Qingbing Zheng, Kiryl D. Piatkevich, Zhifei Fu. A highly stable monomeric red fluorescent protein for advanced microscopy. Nature Methods 2025, 8 https://doi.org/10.1038/s41592-025-02676-5
Helge Vatheuer, Oscar Palomino‐Hernández, Janis Müller, Phillip Galonska, Serghei Glinca, Paul Czodrowski. Protonation Effects in Protein‐Ligand Complexes – A Case Study of Endothiapepsin and Pepstatin A with Computational and Experimental Methods. ChemMedChem 2025, 20 (8) https://doi.org/10.1002/cmdc.202400953
Varun Mandalaparthy, Nico F. A. van der Vegt. A generic model for pH-sensitive collapse of hydrophobic polymers. Physical Chemistry Chemical Physics 2025, 27 (14) , 6984-6993. https://doi.org/10.1039/D4CP04756G
Li-Nan Chen, Hui Zhou, Kun Xi, Shizhuo Cheng, Yongfeng Liu, Yifan Fu, Xiangyu Ma, Ping Xu, Su-Yu Ji, Wei-Wei Wang, Dan-Dan Shen, Huibing Zhang, Qingya Shen, Renjie Chai, Min Zhang, Lin Yang, Feng Han, Chunyou Mao, Xiujun Cai, Yan Zhang. Proton perception and activation of a proton-sensing GPCR. Molecular Cell 2025, 85 (8) , 1640-1657.e8. https://doi.org/10.1016/j.molcel.2025.02.030
Ming Zhang, Minmin Chen, Yizhe Yan, Juan Lu, Jun Sheng, Mingying Gui, Xiao Ma. Comprehensive characterisation of bioactive compounds in Boletus edulis as functional foods to improve muscle atrophy; through whole plant targeted metabolomics, network pharmacology, in vivo and in vitro experiments, molecular docking and molecular dynamics analysis. Journal of Ethnopharmacology 2025, 346 , 119685. https://doi.org/10.1016/j.jep.2025.119685
Luca Manciocchi, Alexandre Bianchi, Valérie Mazan, Mark Potapov, Katharina M. Fromm, Martin Spichty. Modeling the Interaction Between Silver(I) Ion and Proteins with 12-6 Lennard-Jones Potential: A Bottom-Up Parameterization Approach. Biophysica 2025, 5 (1) , 7. https://doi.org/10.3390/biophysica5010007
Vaibhav Modi, Dmitry Morzov. FluProCAD: a computational screening workflow for fluorescent protein variants. Molecular Physics 2025, 2 https://doi.org/10.1080/00268976.2025.2458641
Matthew K. Howard, Nicholas Hoppe, Xi-Ping Huang, Darko Mitrovic, Christian B. Billesbølle, Christian B. Macdonald, Eshan Mehrotra, Patrick Rockefeller Grimes, Donovan D. Trinidad, Lucie Delemotte, Justin G. English, Willow Coyote-Maestas, Aashish Manglik. Molecular basis of proton sensing by G protein-coupled receptors. Cell 2025, 188 (3) , 671-687.e20. https://doi.org/10.1016/j.cell.2024.11.036
Qiuyu Han, Yuxin Chen, Ziyang Zhang, Liping Fan, Jiaoyue Qiu, Wenqing Zhang, Qi Chen, Jinhui Xu, Qianlian Wu, Yue Zhang, Hongbo Liu, Zhishu Tang, Bo Li, Huaxu Zhu. Unveiling the impact of self-assembly on ultrafiltration: Insights from salvianolic acid B. Separation and Purification Technology 2025, 354 , 128857. https://doi.org/10.1016/j.seppur.2024.128857
Ya-Kun Zhang, Zhan Xue, Jian-Bo Tong, Jing Tan, Min Yang, Yan-Rong Zeng, Cheng-Jian Tan. Exploration and computational assessment of ochrocephalamine G from Oxytropis ochrocephala as an anti-HBV candidate. Journal of Asian Natural Products Research 2025, 57 , 1-14. https://doi.org/10.1080/10286020.2024.2441773
Rontu Das, Debashis Kundu. Enlightenment the dynamic behavior of norbornene–modified ’click’ 4–arm polyethylene glycol hydrogel: Delving into framework properties and transport properties through molecular dynamics simulations. Computational Materials Science 2025, 247 , 113516. https://doi.org/10.1016/j.commatsci.2024.113516
Yasmin Momin, Vilas Beloshe. Pharmacophore modeling in drug design. 2025, 313-324. https://doi.org/10.1016/bs.apha.2025.01.010
A. Scacchi, M. Vuorte, M. Sammalkorpi. Multiscale modelling of biopolymers. Advances in Physics: X 2024, 9 (1) https://doi.org/10.1080/23746149.2024.2358196
Ziyuan Niu, Georgios Kementzidis, Peng Zhang, Ziji Zhang, Bernard Essuman, Karin Hasegawa, Miriam Rafailovich, Marcia Simon, Bertal Huseyin Aktas, Evangelos Papadopoulos, Yuefan Deng. Modelling of SARS-CoV-2 spike protein structures at varying pH values. Molecular Simulation 2024, 50 (17-18) , 1540-1552. https://doi.org/10.1080/08927022.2024.2415524
Peter Vogel, David Beyer, Christian Holm, Thomas Palberg. CO 2 -induced drastic decharging of dielectric surfaces in aqueous suspensions. Soft Matter 2024, 20 (46) , 9261-9272. https://doi.org/10.1039/D4SM00957F
Zilan Feng, Chuan Li, Xiangzhou Yi, Changfeng Xue, Xia Gao, Lin Liao, Qiongyao Xiang, Xuanri Shen, Zhisheng Pei. Raman spectroscopy and molecular dynamics simulations of protein microgels at the oil-water interface. International Journal of Biological Macromolecules 2024, 279 , 135398. https://doi.org/10.1016/j.ijbiomac.2024.135398
Hongpeng Yu, Xiaotong Wei, Huan Ding, Shaodan Hu, Feng Sun, Zhenghua Cao, Li Shi. Exploring the potential mechanisms of Mahuang Fuzi Xixin decoction in treating elderly bronchial asthma through network pharmacology, molecular docking, and molecular dynamics simulations. Medicine 2024, 103 (41) , e39921. https://doi.org/10.1097/MD.0000000000039921
Xiaocheng Guo, Xinyuan Feng, Yue Yang, He Zhang, Lunhao Bai. Spermidine attenuates chondrocyte inflammation and cellular pyroptosis through the AhR/NF-κB axis and the NLRP3/caspase-1/GSDMD pathway. Frontiers in Immunology 2024, 15 https://doi.org/10.3389/fimmu.2024.1462777
Dilpreet Singh, Suprit Dilip Saoji. Nanoarchitectonics in Macromolecular Science: Integrating Molecular Dynamics with Smart Materials. Journal of Macromolecular Science, Part B 2024, , 1-10. https://doi.org/10.1080/00222348.2024.2403861
Ya-Kun Zhang, Jian-Bo Tong, Mu-Xuan Luo, Xiao-Yu Xing, Yu-Lu Yang, Zhi-Peng Qing, Ze-Lei Chang, Yan-Rong Zeng. Design and evaluation of piperidine carboxamide derivatives as potent ALK inhibitors through 3D-QSAR modeling, artificial neural network and computational analysis. Arabian Journal of Chemistry 2024, 17 (9) , 105863. https://doi.org/10.1016/j.arabjc.2024.105863
Kelly M. Lee, Vance W. Jaeger. Adsorption of Staphylococcus aureus biofilm associated compounds on silica probed with molecular dynamics simulations. Biointerphases 2024, 19 (5) https://doi.org/10.1116/6.0003870
Alexander Zlobin, Ivan Smirnov, Andrey Golovin. Dynamic interchange between two protonation states is characteristic of active sites of cholinesterases. Protein Science 2024, 33 (8) https://doi.org/10.1002/pro.5100
Simon M Lichtinger, Joanne L Parker, Simon Newstead, Philip C Biggin. The mechanism of mammalian proton-coupled peptide transporters. eLife 2024, 13 https://doi.org/10.7554/eLife.96507
Simon M Lichtinger, Joanne L Parker, Simon Newstead, Philip C Biggin. The mechanism of mammalian proton-coupled peptide transporters. eLife 2024, 13 https://doi.org/10.7554/eLife.96507.3
Carter J. Wilson, Bert L. de Groot, Vytautas Gapsys. Resolving coupled pH titrations using alchemical free energy calculations. Journal of Computational Chemistry 2024, 45 (17) , 1444-1455. https://doi.org/10.1002/jcc.27318
Azadeh Eskandari, Thean Chor Leow, Mohd Basyaruddin Abdul Rahman, Siti Nurbaya Oslan. Molecular dynamics-guided insight into the adsorption–inhibition mechanism for controlling ice growth/melt of antifreeze protein type IV mutant from longhorn sculpin fish. Chemical Papers 2024, 78 (7) , 4437-4454. https://doi.org/10.1007/s11696-024-03407-4
Subramani Selvi, Arun Kannan, John M. Jayaraj, Thangavel Selvi, Muthusamy Karthikeyan, Chidambaram Prahalathan, Natarajan Sampath, Kannupal Srinivasan. Synthesis and biological evaluation of 3‐hydroxypyrazoles as aquaporin 9 inhibitors. Journal of Heterocyclic Chemistry 2024, 61 (4) , 669-679. https://doi.org/10.1002/jhet.4793
Hamed Zahraee, Seyed Shahriar Arab, Zahra Khoshbin, Seyed Mohammad Taghdisi, Khalil Abnous. Molecular dynamics simulation as a promising approach for computational study of liquid crystal-based aptasensors. Journal of Biomolecular Structure and Dynamics 2024, 40 , 1-13. https://doi.org/10.1080/07391102.2024.2315326
Sunita Kumari, Rudolf Podgornik. On the nature of screening in charge-regulated macroion solutions. The Journal of Chemical Physics 2024, 160 (1) https://doi.org/10.1063/5.0187324
Patrick Duchstein, Felix Löffler, Dirk Zahn. Efficient Assessment of ‘Instantaneous pK’ Values from Molecular Dynamics Simulations. ChemPhysChem 2024, 25 (1) https://doi.org/10.1002/cphc.202300489
Oleksandr Khoshaba, Viktor Grechaninov, Tetiana Molodetska, Kostiantyn Zavertailo, Illia Malinich. Structural and Functional Features of the Synthetic Benchmark. 2024, 103-115. https://doi.org/10.1007/978-981-97-7571-2_9
Fernando Luís Barroso da Silva. Constant-pH Simulation Methods for Biomolecular Systems. 2024, 942-963. https://doi.org/10.1016/B978-0-12-821978-2.00090-8
Giovanni Settanni. Computational approaches to lipid-based nucleic acid delivery systems. The European Physical Journal E 2023, 46 (12) https://doi.org/10.1140/epje/s10189-023-00385-5
Hu Ruixuan, Arghya Majee, Jure Dobnikar, Rudolf Podgornik. Electrostatic interactions between charge regulated spherical macroions. The European Physical Journal E 2023, 46 (11) https://doi.org/10.1140/epje/s10189-023-00373-9
Jiaxin Luo, Aoqi Zhang, Yuan Yao, Jun Yuan. Finding small molecular compounds to decrease trimethylamine oxide levels in atherosclerosis by virtual screening. Open Chemistry 2023, 21 (1) https://doi.org/10.1515/chem-2023-0128
Shuai Zhang, Chao Zhang, Min Zhang, Xin Liu, Sheng Xue. Model construction and optimization of coal molecular structure. Journal of Molecular Structure 2023, 1290 , 135960. https://doi.org/10.1016/j.molstruc.2023.135960
Anton O. Chugunov, Elena A. Dvoryakova, Maria A. Dyuzheva, Tatyana R. Simonyan, Valeria F. Tereshchenkova, Irina Yu. Filippova, Roman G. Efremov, Elena N. Elpidina. Fighting Celiac Disease: Improvement of pH Stability of Cathepsin L In Vitro by Computational Design. International Journal of Molecular Sciences 2023, 24 (15) , 12369. https://doi.org/10.3390/ijms241512369
Tuuva Kastinen, Dawid Lupa, Piotr Bonarek, Dmitrii Fedorov, Maria Morga, Markus B. Linder, Jodie L. Lutkenhaus, Piotr Batys, Maria Sammalkorpi. pH dependence of the assembly mechanism and properties of poly( l -lysine) and poly( l -glutamic acid) complexes. Physical Chemistry Chemical Physics 2023, 25 (27) , 18182-18196. https://doi.org/10.1039/D3CP01421E
Lorena Zuzic, Jan K Marzinek, Ganesh S Anand, Jim Warwicker, Peter J Bond. A pH-dependent cluster of charges in a conserved cryptic pocket on flaviviral envelopes. eLife 2023, 12 https://doi.org/10.7554/eLife.82447
Adam L. Harmat, Maria Morga, Jodie L. Lutkenhaus, Piotr Batys, Maria Sammalkorpi. Molecular mechanisms of pH-tunable stability and surface coverage of polypeptide films. Applied Surface Science 2023, 615 , 156331. https://doi.org/10.1016/j.apsusc.2023.156331
Yu Wang, Hao Yang, Wei He, Peixuan Sun, Wenjin Zhao, Miao Liu. Exploring the Potential Hormonal Effects of Tire Polymers (TPs) on Different Species Based on a Theoretical Computational Approach. Polymers 2023, 15 (7) , 1719. https://doi.org/10.3390/polym15071719
Nicolas Pierre Friedrich Mueller, Paolo Carloni, Mercedes Alfonso-Prieto. Molecular determinants of acrylamide neurotoxicity through covalent docking. Frontiers in Pharmacology 2023, 14 https://doi.org/10.3389/fphar.2023.1125871
Pierre Tufféry, Philippe Derreumaux. A refined pH-dependent coarse-grained model for peptide structure prediction in aqueous solution. Frontiers in Bioinformatics 2023, 3 https://doi.org/10.3389/fbinf.2023.1113928
Fang-Fang Luo, Zhi-Tao Cai, Yan-Dong Huang, . Progress in protein p<i>K</i><sub>a</sub> prediction. Acta Physica Sinica 2023, 72 (24) , 248704. https://doi.org/10.7498/aps.72.20231356
Luca Donati, Marcus Weber. Assessing transition rates as functions of environmental variables. The Journal of Chemical Physics 2022, 157 (22) https://doi.org/10.1063/5.0109555
Lisbeth R. Kjølbye, Gilberto P. Pereira, Alessio Bartocci, Martina Pannuzzo, Simone Albani, Alessandro Marchetto, Brian Jiménez-García, Juliette Martin, Giulia Rossetti, Marco Cecchini, Sangwook Wu, Luca Monticelli, Paulo C. T. Souza. Towards design of drugs and delivery systems with the Martini coarse-grained model. QRB Discovery 2022, 3 https://doi.org/10.1017/qrd.2022.16

Download PDF

Get e-Alerts

Journal of Chemical Theory and Computation

Cite this: J. Chem. Theory Comput. 2022, 18, 10, 6148–6160

Click to copy citationCitation copied!

https://doi.org/10.1021/acs.jctc.2c00516

Published September 21, 2022

CC-BY 4.0 .

Article Views

Altmetric

Citations

Learn about these metrics

Article Views are the COUNTER-compliant sum of full text article downloads since November 2008 (both PDF and HTML) across all institutions and individuals. These metrics are regularly updated to reflect usage leading up to the last few days.

Citations are the number of other articles citing this article, calculated by Crossref and updated daily. Find more information about Crossref citation counts.

The Altmetric Attention Score is a quantitative measure of the attention that a research article has received online. Clicking on the donut icon will load a page at altmetric.com with additional details about the score and the social media presence for the given article. Find more information on the Altmetric Attention Score and how the score is calculated.

Abstract
High Resolution Image
Download MS PowerPoint Slide
Figure 1
Figure 1. Illustration of the additional λ-dependent potential energy terms (B–D). Panel A shows the protonation free energy of a titratable residue in its reference state obtained at the force field level, ΔG_i^MM(λ_i) . To compensate for the shortcomings of the force field and obtain a zero free energy difference between the two protonation states (A + B), we add the correction potential, V_i^MM(λ), shown in panel B. A biasing potential, V_i^bias(λ), (42) is introduced to avoid sampling of nonphysical states (panel C). To model the proton chemical potential (pH), we add a pH-dependent term, V^pH(λ_i) (panel D). For a titratable residue at pH ≠ pK_a, the total λ-dependent potential, including the interpolated force field functions and the three additional terms, is illustrated in panel (A + B + C + D).
High Resolution Image
Download MS PowerPoint Slide
Figure 2
Figure 2. Multisite representation illustrated for a histidine side chain. Each possible protonation state is represented by its own λ-coordinate. HSP refers to doubly protonated histidine, HSD, and HSE to histidine singly protonated at the N_δ or N_ϵ, respectively. HS0 is a common, nonphysical state of the residue. To restrict the sampling to the plane connecting the physical states, a constraint λ₁ + λ₂ + λ₃ = 1 is applied (gray plane). A biasing potential is also applied to enhance sampling at the end states, where one of the λ-coordinates is close to one, while the other coordinates are close to zero.
High Resolution Image
Download MS PowerPoint Slide
Figure 3
Figure 3. Titrations of tripeptide amino acids (Glu, Asp, and His) in water. The top and bottom rows show titrations performed with the modified AA CHARMM36 (55) and CG Martini (41) force fields, respectively. In all simulations, neutrality was maintained by including ten buffer particles in combination with the charge constraint. Dots show the fraction of frames in which the residue is deprotonated, and the dashed lines represent the fits to the Henderson–Hasselbalch equation. For His, the blue color represents the macroscopic pK_a, while yellow and red represent the microscopic pK_a values for HSD (proton on N_δ) and HSE (proton on N_ϵ), respectively. In the Martini 2.0 model, HSD and HSE are indistinguishable and hence only the macroscopic titration curve is shown. Errors of S^deprot were estimated from the standard error of the mean for the different replicas. From the fits, the pK_a values were estimated and listed in Table 1.
High Resolution Image
Download MS PowerPoint Slide
Figure 4
Figure 4. Titration curves of the cardiotoxin V protein obtained from constant pH MD simulations with the CHARMM36 (top) and Martini 2.0 force fields (bottom). For each of the four titratable residues in this protein, the dots show the fraction of frames in which the residue is deprotonated. Errors of S^deprot were estimated from the standard error of the mean for the different replicas. The lines show the best fits to the Henderson–Hasselbalch equation. The pK_a values for each titratable residue were estimated from these fits and listed in Table 1. The right panel shows the protein structure with the four titratable residues highlighted in stick representation.
High Resolution Image
Download MS PowerPoint Slide
Figure 5
Figure 5. Titration curves of the HEWL protein obtained from constant pH MD simulations with the CHARMM36 (top) and Martini 2.0 force fields (bottom). For each of the ten titratable residues, the dots show the fraction of frames in which that residue is deprotonated. Errors of S^deprot were estimated from the standard error of the mean for the different replicas. The lines show the best fit to the Henderson-Hasselbach equation. The pK_a values for each titratable residue were estimated from these fits and listed in Table 1. The right panel shows the protein structure with the ten titratable residues highlighted in stick representation.
High Resolution Image
Download MS PowerPoint Slide
Figure 6
Figure 6. (A) Relative performance of interpolating potentials in a previous implementation of CpHMD into a fork of GROMACS 3.3 release (blue) and of charge interpolation in our new implementation (red) as a function of the number of titratable sites. The simulations were performed for the turkey ovomucoid inhibitor protein (PDB ID: 2GKR (53)), shown in the inset, where the titratable sites are highlighted in stick representation. (B) Comparison of the performance between CPU-only and CPU+GPU implementations for the ligand-gated ion channel GLIC (PDB ID: 4HFI (52)) with 185 titratable sites. In total, the GLIC system contained 292135 atoms.
High Resolution Image
Download MS PowerPoint Slide
References
This article references 76 other publications.
1. 1
  Hollingsworth, S. A.; Dror, R. O. Molecular Dynamics Simulation for All. Neuron 2018, 99, 1129– 1143, DOI: 10.1016/j.neuron.2018.08.011
  1
  Molecular Dynamics Simulation for All
  Hollingsworth, Scott A.; Dror, Ron O.
  Neuron (2018), 99 (6), 1129-1143CODEN: NERNET; ISSN:0896-6273. (Cell Press)
  A review. The impact of mol. dynamics (MD) simulations in mol. biol. and drug discovery has expanded dramatically in recent years. These simulations capture the behavior of proteins and other biomols. in full at. detail and at very fine temporal resoln. Major improvements in simulation speed, accuracy, and accessibility, together with the proliferation of exptl. structural data, have increased the appeal of biomol. simulation to experimentalists-a trend particularly noticeable in, although certainly not limited to, neuroscience. Simulations have proven valuable in deciphering functional mechanisms of proteins and other biomols., in uncovering the structural basis for disease, and in the design and optimization of small mols., peptides, and proteins. Here we describe, in practical terms, the types of information MD simulations can provide and the ways in which they typically motivate further exptl. work.
2. 2
  Groenhof, G.; Modi, V.; Morozov, D. Observe while it happens: catching photoactive proteins in the act with non-adiabatic molecular dynamics simulations. Curr. Opin. Struct. Biol. 2020, 61, 106– 112, DOI: 10.1016/j.sbi.2019.12.013
  2
  Observe while it happens: catching photoactive proteins in the act with non-adiabatic molecular dynamics simulations
  Groenhof, Gerrit; Modi, Vaibhav; Morozov, Dmitry
  Current Opinion in Structural Biology (2020), 61 (), 106-112CODEN: COSBEF; ISSN:0959-440X. (Elsevier Ltd.)
  A review. Organisms use photo-receptors to react to light. The first step is usually the absorption of a photon by a prosthetic group embedded inside the photo-receptor, often a conjugated chromophore. The electronic changes in the chromophore induced by photo-absorption can trigger a cascade of structural or chem. transformations that culminate into a response to light. Understanding how these proteins have evolved to mediate their activation process has remained challenging because the required time and spacial resolns. are notoriously difficult to achieve exptl. Therefore, mechanistic insights into photoreceptor activation have been predominantly obtained with computer simulations. Here we briefly outline the challenges assocd. with such computations and review the progress made in this field.
3. 3
  Warshel, A.; Sussman, F.; King, G. Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process. Biochemistry 1986, 25, 8368– 8372, DOI: 10.1021/bi00374a006
  3
  Free energy of charges in solvated proteins: microscopic calculations using a reversible charging process
  Warshel, A.; Sussman, F.; King, G.
  Biochemistry (1986), 25 (26), 8368-72CODEN: BICHAW; ISSN:0006-2960.
  Evaluation of the free energy of ionization of acidic groups in proteins may be used as a powerful and general test case for detg. the reliability of calcns. of electrostatic energies in macromols. This test case was addressed by using an adiabatic charging process that evaluates the changes in free energies assocd. with ionizing the acidic groups of aspartate-3 and glutamate-7 in bovine pancreatic trypsin inhibitor and aspartic acid in soln. The results of these free energy calcns. are very encouraging; the error range is ∼1 kcal/mol for these values, which are ∼70 kcal/mol. Thus, the stage of obtaining quant. results in modeling the energetics of solvated proteins is finally being approached.
4. 4
  Alexov, E.; Mehler, E. L.; Baker, N.; Baptista, A. M.; 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 pKa values in proteins. Proteins: Struct., Funct., Bioinf. 2011, 79, 3260– 3275, DOI: 10.1002/prot.23189
  4
  Progress in the prediction of pKa values in proteins
  Alexov, Emil; Mehler, Ernest L.; Baker, Nathan; Baptista, Antonio M.; Huang, Yong; Milletti, Francesca; Erik Nielsen, Jens; Farrell, Damien; Carstensen, Tommy; Olsson, Mats H. M.; Shen, Jana K.; Warwicker, Jim; Williams, Sarah; Word, J. Michael
  Proteins: Structure, Function, and Bioinformatics (2011), 79 (12), 3260-3275CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
  A review. The pKa-cooperative aims to provide a forum for exptl. and theor. researchers interested in protein pKa values and protein electrostatics in general. The first round of the pKa-cooperative, which challenged computational labs to carry out blind predictions against pKas exptl. detd. in the lab. of Bertrand Garcia-Moreno, was completed and results discussed at the Telluride meeting (July 6-10, 2009). This article serves as an introduction to the reports submitted by the blind prediction participants that will be published in a special issue of PROTEINS: Structure, Function and Bioinformatics. Here, the authors briefly outline existing approaches for pKa calcns., emphasizing methods that were used by the participants in calcg. the blind pKa values in the first round of the cooperative. The authors then point out some of the difficulties encountered by the participating groups in making their blind predictions, and finally try to provide some insights for future developments aimed at improving the accuracy of pKa calcns. Proteins 2011; © 2011 Wiley-Liss, Inc.
5. 5
  Chen, W.; Morrow, B. H.; Shi, C.; Shen, J. K. Recent development and application of constant pH molecular dynamics. Mol. Simul. 2014, 40, 830– 838, DOI: 10.1080/08927022.2014.907492
  5
  Recent development and application of constant pH molecular dynamics
  Chen, Wei; Morrow, Brian H.; Shi, Chuanyin; Shen, Jana K.
  Molecular Simulation (2014), 40 (10-11), 830-838CODEN: MOSIEA; ISSN:0892-7022. (Taylor & Francis Ltd.)
  A review. Soln. pH is a crit. environmental factor for chem. and biol. processes. Over the last decade, significant efforts have been made in the development of const. pH mol. dynamics (pHMD) techniques for gaining detailed insights into pH-coupled dynamical phenomena. In this article, we review the advancement of this field in the past five years, placing a special emphasis on the development of the all-atom continuous pHMD technique. We discuss various applications, including the prediction of pKa shifts for proteins, nucleic acids and surfactant assemblies, elucidation of pH-dependent population shifts, protein-protein and protein-RNA binding, as well as the mechanisms of pH-dependent self-assembly and phase transitions of surfactants and peptides. We also discuss future directions for the further improvement of the pHMD techniques.
6. 6
  Zeng, X.; Mukhopadhyay, S.; Brooks, C. L., III Residue-level resolution of alphavirus envelope protein interactions in pH-dependent fusion. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 2034– 2039, DOI: 10.1073/pnas.1414190112
  6
  Residue-level resolution of alphavirus envelope protein interactions in pH-dependent fusion
  Zeng, Xiancheng; Mukhopadhyay, Suchetana; Brooks, Charles L., III
  Proceedings of the National Academy of Sciences of the United States of America (2015), 112 (7), 2034-2039CODEN: PNASA6; ISSN:0027-8424. (National Academy of Sciences)
  Alphavirus envelope proteins, organized as trimers of E2-E1 heterodimers on the surface of the pathogenic alphavirus, mediate the low pH-triggered fusion of viral and endosomal membranes in human cells. The lack of specific treatment for alphaviral infections motivates our exploration of potential antiviral approaches by inhibiting one or more fusion steps in the common endocytic viral entry pathway. In this work, we performed const. pH mol. dynamics based on an at. model of the alphavirus envelope with icosahedral symmetry. We have identified pH-sensitive residues that cause the largest shifts in thermodn. driving forces under neutral and acidic pH conditions for various fusion steps. A series of conserved interdomain His residues is identified to be responsible for the pH-dependent conformational changes in the fusion process, and ligand binding sites in their vicinity are anticipated to be potential drug targets aimed at inhibiting viral infections.
7. 7
  Law, S. M.; Zhang, B. W.; Brooks, C. L., III pH-sensitive residues in the p19 RNA silencing suppressor protein from carnation Italian ringspot virus affect siRNA binding stability. Protein Sci. 2013, 22, 595– 604, DOI: 10.1002/pro.2243
  7
  pH-sensitive residues in the p19 RNA silencing suppressor protein from carnation Italian ringspot virus affect siRNA binding stability
  Law, Sean M.; Zhang, Bin W.; Brooks, Charles L., III
  Protein Science (2013), 22 (5), 595-604CODEN: PRCIEI; ISSN:1469-896X. (Wiley-Blackwell)
  Tombusviruses, such as Carnation Italian ringspot virus (CIRV), encode a protein homodimer called p19 that is capable of suppressing RNA silencing in their infected hosts by binding to and sequestering short-interfering RNA (siRNA) away from the RNA silencing pathway. P19 binding stability has been shown to be sensitive to changes in pH but the specific amino acid residues involved have remained unclear. Using const. pH mol. dynamics simulations, we have identified key pH-dependent residues that affect CIRV p19-siRNA binding stability at various pH ranges based on calcd. changes in the free energy contribution from each titratable residue. At high pH, the deprotonation of Lys60, Lys67, Lys71, and Cys134 has the largest effect on the binding stability. Similarly, deprotonation of several acidic residues (Asp9, Glu12, Asp20, Glu35, and/or Glu41) at low pH results in a decrease in binding stability. At neutral pH, residues Glu17 and His132 provide a small increase in the binding stability and we find that the optimal pH range for siRNA binding is between 7.0 and 10.0. Overall, our findings further inform recent expts. and are in excellent agreement with data on the pH-dependent binding profile.
8. 8
  Ellis, C. R.; Shen, J. pH-dependent population shift regulates BACE1 activity and inhibition. J. Am. Chem. Soc. 2015, 137, 9543– 9546, DOI: 10.1021/jacs.5b05891
  8
  pH-Dependent Population Shift Regulates BACE1 Activity and Inhibition
  Ellis, Christopher R.; Shen, Jana
  Journal of the American Chemical Society (2015), 137 (30), 9543-9546CODEN: JACSAT; ISSN:0002-7863. (American Chemical Society)
  BACE1, a major therapeutic target for treatment of Alzheimer's disease, functions within a narrow pH range. Despite tremendous effort and progress in the development of BACE1 inhibitors, details of the underlying pH-dependent regulatory mechanism remain unclear. Here the authors elucidate the pH-dependent conformational mechanism that regulates BACE1 activity using continuous const.-pH mol. dynamics (MD). The simulations reveal that BACE1 mainly occupies three conformational states and that the relative populations of the states shift according to pH. At intermediate pH, when the catalytic dyad is monoprotonated, a binding-competent state is highly populated, while at low and high pH a Tyr-inhibited state is dominant. The authors' data provide strong evidence supporting conformational selection as a major mechanism for substrate and peptide-inhibitor binding. These new insights, while consistent with expt., greatly extend the knowledge of BACE1 and have implications for further optimization of inhibitors and understanding potential side effects of targeting BACE1. Finally, the work highlights the importance of properly modeling protonation states in MD simulations.
9. 9
  Kim, M. O.; Blachly, P. G.; McCammon, J. A. Conformational dynamics and binding free energies of inhibitors of BACE-1: from the perspective of protonation equilibria. PLoS computational biology 2015, 11, e1004341, DOI: 10.1371/journal.pcbi.1004341
  9
  Conformational dynamics and binding free energies of inhibitors of BACE-1: from the perspective of protonation equilibria
  Kim, M. Olivia; Blachly, Patrick G.; McCammon, J. Andrew
  PLoS Computational Biology (2015), 11 (10), e1004341/1-e1004341/28CODEN: PCBLBG; ISSN:1553-7358. (Public Library of Science)
  BACE-1 is the β-secretase responsible for the initial amyloidogenesis in Alzheimer's disease, catalyzing hydrolytic cleavage of substrate in a pH-sensitive manner. The catalytic mechanism of BACE-1 requires water-mediated proton transfer from aspartyl dyad to the substrate, as well as structural flexibility in the flap region. Thus, the coupling of protonation and conformational equil. is essential to a full in silico characterization of BACE-1. In this work, we perform const. pH replica exchange mol. dynamics simulations on both apo BACE-1 and five BACE-1-inhibitor complexes to examine the effect of pH on dynamics and inhibitor binding properties of BACE-1. In our simulations, we find that soln. pH controls the conformational flexibility of apo BACE-1, whereas bound inhibitors largely limit the motions of the holo enzyme at all levels of pH. The microscopic pKa values of titratable residues in BACE-1 including its aspartyl dyad are computed and compared between apo and inhibitor-bound states. Changes in protonation between the apo and holo forms suggest a thermodn. linkage between binding of inhibitors and protons localized at the dyad. Utilizing our recently developed computational protocol applying the binding polynomial formalism to the const. pH mol. dynamics (CpHMD) framework, we are able to obtain the pH-dependent binding free energy profiles for various BACE-1-inhibitor complexes. Our results highlight the importance of correctly addressing the binding-induced protonation changes in protein-ligand systems where binding accompanies a net proton transfer. This work comprises the first application of our CpHMD-based free energy computational method to protein-ligand complexes and illustrates the value of CpHMD as an all-purpose tool for obtaining pH-dependent dynamics and binding free energies of biol. systems.
10. 10
  Sarkar, A.; Gupta, P. L.; Roitberg, A. E. pH-dependent conformational changes due to ionizable residues in a hydrophobic protein interior: The study of L25K and L125K variants of SNase. J. Phys. Chem. B 2019, 123, 5742– 5754, DOI: 10.1021/acs.jpcb.9b03816
  10
  pH-Dependent Conformational Changes Due to Ionizable Residues in a Hydrophobic Protein Interior: The Study of L25K and L125K Variants of SNase
  Sarkar, Ankita; Gupta, Pancham Lal; Roitberg, Adrian E.
  Journal of Physical Chemistry B (2019), 123 (27), 5742-5754CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
  Ionizable residues in the hydrophobic interior of certain proteins are known to play important roles in life processes like energy transduction and enzyme catalysis. These internal ionizable residues show exptl. apparent pKa values having large shifts as compared to their values in soln. In the present work, we study the pH-dependent conformational changes undergone by two variants of staphylococcal nuclease (SNase), L25K and L125K, using pH replica exchange mol. dynamics (pH-REMD) in explicit solvent. Our results show that the obsd. pKa of Lys25 and Lys125 are significantly different than their pKa in soln. We obsd. that the internal lysine residues prefer to be water-exposed when protonated at low pH, but they remain buried within the hydrophobic pocket when deprotonated at high pH. Using thermodn. laws, we est. the microscopic conformation-specific pKa of the water-exposed and buried conformations of the internal lysine residues and explain their relation to the macroscopic obsd. pKa values. We present the differences in the microscopic mechanisms that lead to similar exptl. obsd. apparent pKa of Lys25 and Lys125, and explain the need of thermodn. models of different complexities to account for our calcns. We see that L25K displays pH-dependent fluctuations throughout the entire β barrel and the α1 helix. In contrast, pH-independent fluctuations are obsd. in L125K, primarily limited to the α3 helix. The present computational study offers a detailed atomistic understanding of the determinants of the obsd. anomalous pKa of internal ionizable residues, bolstering the exptl. findings.
11. 11
  Sarkar, A.; Roitberg, A. E. PH-Dependent Conformational Changes Lead to a Highly Shifted p K a for a Buried Glutamic Acid Mutant of SNase. J. Phys. Chem. B 2020, 124, 11072– 11080, DOI: 10.1021/acs.jpcb.0c07136
  11
  pH-Dependent Conformational Changes Lead to a Highly Shifted pKa for a Buried Glutamic Acid Mutant of SNase
  Sarkar, Ankita; Roitberg, Adrian E.
  Journal of Physical Chemistry B (2020), 124 (49), 11072-11080CODEN: JPCBFK; ISSN:1520-5207. (American Chemical Society)
  Ionizable residues are rarely present in the hydrophobic interior of proteins, but when they are, they play important roles in biol. processes such as energy transduction and enzyme catalysis. Internal ionizable residues have anomalous exptl. pKa values with respect to their pKa in bulk water. This work studies the atomistic cause of the highly shifted pKa of the internal Glu23 in the artificially mutated variant V23E of Staphylococcal Nuclease (SNase) using pH replica exchange mol. dynamics (pH-REMD) simulations. The pKa of Glu23 obtained from the authors' calcns. is 6.55, which is elevated with respect to the glutamate pKa of 4.40 in bulk water. The calcd. value is close to the exptl. pKa of 7.10. The authors' simulations show that the highly shifted pKa of Glu23 is the product of a pH-dependent conformational change, which was obsd. exptl. and also seen in the authors' simulations. The authors carry out an anal. of this pH-dependent conformational change in response to the protonation state change of Glu23. Using a four-state thermodn. model, the authors est. the two conformation-specific pKa values of Glu23 and describe the coupling between the conformational and ionization equil.
12. 12
  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: Struct., Funct., Bioinf. 1997, 27, 523– 544, DOI: 10.1002/(SICI)1097-0134(199704)27:4<523::AID-PROT6>3.0.CO;2-B
  12
  Simulation of protein conformational freedom as a function of pH: constant-pH molecular dynamics using implicit titration
  Baptista, Antonio M.; Martel, Paulo J.; Petersen, Steffen B.
  Proteins: Structure, Function, and Genetics (1997), 27 (4), 523-544CODEN: PSFGEY; ISSN:0887-3585. (Wiley-Liss)
  Soln. pH is a determinant parameter on protein function and stability, and its inclusion in mol. dynamics simulations is attractive for studies at the mol. level. Current mol. dynamics simulations can consider pH only in a very limited way, through a somewhat arbitrary choice of a set of fixed charges on the titrable sites. Conversely, continuum electrostatic methods that explicitly treat pH effects assume a single protein conformation whose choice is not clearly defined. In this paper we describe a general method that combines both titrn. and conformational freedom. The method is based on a potential of mean force for implicit titrn. and combines both usual mol. dynamics and pH-dependent calcns. based on continuum methods. A simple implementation of the method, using a mean field approxn., is presented and applied to the bovine pancreatic trypsin inhibitor. We believe that this const.-pH mol. dynamics method, by correctly sampling both charges and conformation, can become a valuable help in the understanding of the dependence of protein function and stability on pH.
13. 13
  Baptista, A. M.; Teixeira, V. H.; Soares, C. M. Constant-pH molecular dynamics using stochastic titration. J. Chem. Phys. 2002, 117, 4184– 4200, DOI: 10.1063/1.1497164
  13
  Constant-pH molecular dynamics using stochastic titration
  Baptista, Antonio M.; Teixeira, Vitor H.; Soares, Claudio M.
  Journal of Chemical Physics (2002), 117 (9), 4184-4200CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  A new method is proposed for performing const.-pH mol. dynamics (MD) simulations, i.e., MD simulations where pH is one of the external thermodn. parameters, like the temp. or the pressure. The protonation state of each titrable site in the solute is allowed to change during a mol. mechanics (MM) MD simulation, the new states being obtained from a combination of continuum electrostatics (CE) calcns. and Monte Carlo (MC) simulation of protonation equil. The coupling between the MM/MD and CE/MC algorithms is done in a way that ensures a proper Markov chain, sampling from the intended semigrand canonical distribution. This stochastic titrn. method is applied to succinic acid, aimed at illustrating the method and examg. the choice of its adjustable parameters. The complete titrn. of succinic acid, using const.-pH MD simulations at different pH values, gives a clear picture of the coupling between the trans/gauche isomerization and the protonation process, making it possible to reconcile some apparently contradictory results of previous studies. The present const.-pH MD method is shown to require a moderate increase of computational cost when compared to the usual MD method.
14. 14
  Bürgi, R.; Kollman, P. A.; van Gunsteren, W. F. Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation. Proteins: Struct., Funct., Bioinf. 2002, 47, 469– 480, DOI: 10.1002/prot.10046
  14
  Simulating proteins at constant pH: an approach combining molecular dynamics and Monte Carlo simulation
  Burgi, Roland; Kollman, Peter A.; Van Gunsteren, Wilfred F.
  Proteins: Structure, Function, and Genetics (2002), 47 (4), 469-480CODEN: PSFGEY; ISSN:0887-3585. (Wiley-Liss, Inc.)
  For the structure and function of proteins, the pH of the soln. is one of the detg. parameters. Current mol. dynamics (MD) simulations account for the soln. pH only in a limited way by keeping each titratable site in a chosen protonation state. We present an algorithm that generates trajectories at a Boltzmann distributed ensemble of protonation states by a combination of MD and Monte Carlo (MC) simulation. The algorithm is useful for pH-dependent structural studies and to investigate in detail the titrn. behavior of proteins. The method is tested on the acidic residues of the protein hen egg white lysozyme. It is shown that small structural changes may have a big effect on the pKA values of titratable residues.
15. 15
  Mongan, J.; Case, D. A.; McCammon, J. A. Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem. 2004, 25, 2038– 2048, DOI: 10.1002/jcc.20139
  15
  Constant pH molecular dynamics in generalized Born implicit solvent
  Mongan, John; Case, David A.; McCammon, J. Andrew
  Journal of Computational Chemistry (2004), 25 (16), 2038-2048CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
  A new method is proposed for const. pH mol. dynamics (MD), employing generalized Born (GB) electrostatics. Protonation states are modeled with different charge sets, and titrating residues sample a Boltzmann distribution of protonation states as the simulation progresses, using Monte Carlo sampling based on GB-derived energies. The method is applied to four different crystal structures of hen egg-white lysozyme (HEWL). PKa predictions derived from the simulations have root-mean-square (RMS) error of 0.82 relative to exptl. values. Similarity of results between the four crystal structures shows the method to be independent of starting crystal structure; this is in contrast to most electrostatics-only models. A strong correlation between conformation and protonation state is noted and quant. analyzed, emphasizing the importance of sampling protonation states in conjunction with dynamics.
16. 16
  Meng, Y.; Roitberg, A. E. Constant pH Replica Exchange Molecular Dynamics in Biomolecules Using a Discrete Protonation Model. J. Chem. Theory Comput. 2010, 6, 1401– 1412, DOI: 10.1021/ct900676b
  16
  Constant pH Replica Exchange Molecular Dynamics in Biomolecules Using a Discrete Protonation Model
  Meng, Yilin; Roitberg, Adrian E.
  Journal of Chemical Theory and Computation (2010), 6 (4), 1401-1412CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  A const. pH replica exchange mol. dynamics (REMD) method is proposed and implemented to improve coupled protonation and conformational state sampling. By mixing conformational sampling at const. pH (with discrete protonation states) with a temp. ladder, this method avoids conformational trapping. Our method was tested and applied to seven different biol. systems. The const. pH REMD not only predicted pKa correctly for small, model compds. but also converged faster than const. pH mol. dynamics (MD). We further tested our const. pH REMD on a heptapeptide from the ovomucoid third domain (OMTKY3). Although const. pH REMD and MD produced very close pKa values, the const. pH REMD showed its advantage in the efficiency of conformational and protonation state samplings.
17. 17
  Itoh, S. G.; Damjanovic, A.; Brooks, B. R. pH replica-exchange method based on discrete protonation states. Proteins: Struct., Funct., Bioinf. 2011, 79, 3420– 3436, DOI: 10.1002/prot.23176
  17
  pH replica-exchange method based on discrete protonation states
  Itoh, Satoru G.; Damjanovic, Ana; Brooks, Bernard R.
  Proteins: Structure, Function, and Bioinformatics (2011), 79 (12), 3420-3436CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
  The authors propose a new algorithm for obtaining proton titrn. curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the const. pH algorithm of J. Mongan et al. (2004). In the original replica-exchange method, simulations of different replicas are performed at different temps., and the temps. are exchanged between the replicas. In the authors' PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original const. pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coeffs. of titratable residues. By performing multiple sets of const. pH simulations started with different initial states, the accuracy of predicted pKa values and Hill coeffs. obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titrn. points and thus better precision of measured pKa values and Hill coeffs. In addn., PHREM exhibits faster convergence of individual simulations than the original const. pH algorithm. Proteins 2011; © 2011 Wiley-Liss, Inc.
18. 18
  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, DOI: 10.1021/ct401042b
  18
  Constant pH Replica Exchange Molecular Dynamics in Explicit Solvent Using Discrete Protonation States: Implementation, Testing, and Validation
  Swails, Jason M.; York, Darrin M.; Roitberg, Adrian E.
  Journal of Chemical Theory and Computation (2014), 10 (3), 1341-1352CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  By utilizing Graphics Processing Units, we show that const. pH mol. dynamics simulations (CpHMD) run in Generalized Born (GB) implicit solvent for long time scales can yield poor pKa predictions as a result of sampling unrealistic conformations. To address this shortcoming, we present a method for performing const. pH mol. dynamics simulations (CpHMD) in explicit solvent using a discrete protonation state model. The method involves std. mol. dynamics (MD) being propagated in explicit solvent followed by protonation state changes being attempted in GB implicit solvent at fixed intervals. Replica exchange along the pH-dimension (pH-REMD) helps to obtain acceptable titrn. behavior with the proposed method. We analyzed the effects of various parameters and settings on the titrn. behavior of CpHMD and pH-REMD in explicit solvent, including the size of the simulation unit cell and the length of the relaxation dynamics following protonation state changes. We tested the method with the amino acid model compds., a small pentapeptide with two titratable sites, and hen egg white lysozyme (HEWL). The proposed method yields superior predicted pKa values for HEWL over hundreds of nanoseconds of simulation relative to corresponding predicted values from simulations run in implicit solvent.
19. 19
  Swails, J. M.; Roitberg, A. E. Enhancing conformation and protonation state sampling of hen egg white lysozyme using pH replica exchange molecular dynamics. J. Chem. Theory Comput. 2012, 8, 4393– 4404, DOI: 10.1021/ct300512h
  19
  Enhancing Conformation and Protonation State Sampling of Hen Egg White Lysozyme Using pH Replica Exchange Molecular Dynamics
  Swails, Jason M.; Roitberg, Adrian E.
  Journal of Chemical Theory and Computation (2012), 8 (11), 4393-4404CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  We evaluate the efficiency of the pH replica exchange mol. dynamics (pH-REMD) method proposed by Itoh et al. (Proteins 2011, 79, 3420-3436) by using it to predict the pKa values of the titratable residues in hen egg white lysozyme (HEWL). PKa values predicted using pH-REMD converge significantly faster than those calcd. using const. pH mol. dynamics (CpHMD). Furthermore, increasing the frequency between exchange attempts in pH-REMD simulations improves protonation and conformational state sampling. By enabling the simulation to sample both conformational and protonation states more rapidly, pH-REMD simulations provide valuable insight into the pH-dependence of HEWL that the CpHMD simulations failed to capture. We present an efficient and highly scalable implementation of pH-REMD as an attractive enhancement to traditional CpHMD methods.
20. 20
  Börjesson, U.; Hünenberger, P. H. Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines. J. Chem. Phys. 2001, 114, 9706– 9719, DOI: 10.1063/1.1370959
  20
  Explicit-solvent molecular dynamics simulation at constant pH: Methodology and application to small amines
  Borjesson, Ulf; Hunenberger, Philippe H.
  Journal of Chemical Physics (2001), 114 (22), 9706-9719CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  A method is developed for performing classical explicit-solvent mol. dynamics (MD) simulations at const. pH, where the protonation state of each ionizable (titratable) group in a simulated compd. is allowed to fluctuate in time, depending on the instantaneous system configuration and the imposed pH. In this method, each ionizable group is treated as a mixed state, i.e., the interaction-function parameters for the group are a linear combination of those of the protonated state and those of the deprotonated state. Free protons are not handled explicitly. Instead, the extent of deprotonation of each group is relaxed towards its equil. value by weak coupling to a "proton bath. " The method relies on precalibrated empirical functions, one for each type of ionizable group present in the simulated compd., which are obtained through multiple MD simulations of monofunctional model compds. In this study, the method is described in detail and its application illustrated by a series of const.-pH MD simulations of small monofunctional amines. In particular, we investigate the influence of the relaxation time used in the weak-coupling scheme, the choice of appropriate model compds. for the calibration of the required empirical functions, and corrections for finite-size effects linked with the small size of the simulation box.
21. 21
  Lee, M. S.; Salsbury, F. R., Jr.; Brooks, C. L., III Constant-pH molecular dynamics using continuous titration coordinates. Proteins: Struct., Funct., Bioinf. 2004, 56, 738– 752, DOI: 10.1002/prot.20128
  21
  Constant-pH molecular dynamics using continuous titration coordinates
  Lee, Michael S.; Salsbury, Freddie R., Jr.; Brooks, Charles L., III
  Proteins: Structure, Function, and Bioinformatics (2004), 56 (4), 738-752CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
  In this work, we explore the question of whether pKa calcns. based on a microscopic description of the protein and a macroscopic description of the solvent can be implemented to examine conformationally dependent proton shifts in proteins. To this end, we introduce a new method for performing const.-pH mol. dynamics (PHMD) simulations utilizing the generalized Born implicit solvent model. This approach employs an extended Hamiltonian in which continuous titrn. coordinates propagate simultaneously with the at. motions of the system. The values adopted by these coordinates are modulated by potentials of mean force of isolated titratable model groups and the pH to control the proton occupation at particular sites in the polypeptide. Our results for four different proteins yield an abs. av. error of ∼1.6 pK units, and point to the role that thermally driven relaxation of the protein environment in the vicinity of titrating groups plays in modulating the local pKa, thereby influencing the obsd. pK1/2 values. While the accuracy of our method is not yet equiv. to methods that obtain pK1/2 values through the ad hoc scaling of electrostatics, the present approach and const. pH methods in general provide a useful framework for studying pH-dependent phenomena. Further work to improve our model to approach quant. agreement with expt. is outlined.
22. 22
  Khandogin, J.; Brooks, C. L. Constant pH Molecular Dynamics with proton tautomerism. Biophys. J. 2005, 89, 141– 157, DOI: 10.1529/biophysj.105.061341
  22
  Constant pH molecular dynamics with proton tautomerism
  Khandogin, Jana; Brooks, Charles L., III
  Biophysical Journal (2005), 89 (1), 141-157CODEN: BIOJAU; ISSN:0006-3495. (Biophysical Society)
  The current article describes a new two-dimensional λ-dynamics method to include proton tautomerism in continuous const. pH mol. dynamics (CPHMD) simulations. The two-dimensional λ-dynamics framework is used to devise a tautomeric state titrn. model for the CPHMD simulations involving carboxyl and histidine residues. Combined with the GBSW implicit solvent model, the new method is tested on titrn. simulations of blocked histidine and aspartic acid as well as two benchmark proteins, turkey ovomucoid third domain (OMTKY3) and RNase A (RNase A). A detailed anal. of the errors inherent to the CPHMD methodol. as well as those due to the underlying solvation model is given. The av. abs. error for the computed pKa values in OMTKY3 is 1.0 pK unit. In RNase A the av. abs. errors for the carboxyl and histidine residues are 1.6 and 0.6 pK units, resp. In contrast to the previous work, the new model predicts the correct sign for all the pKa shifts, but one, in the benchmark proteins. The predictions of the tautomeric states of His12 and His48 and the conformational states of His48 and His119 are in agreement with expt. Based on the simulations of OMTKY3 and RNase A, the current work has demonstrated the capability of the CPHMD technique in revealing pH-coupled conformational dynamics of protein side chains.
23. 23
  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, DOI: 10.1021/ct200061r
  23
  Constant pH Molecular Dynamics in Explicit Solvent with λ-Dynamics
  Donnini, Serena; Tegeler, Florian; Groenhof, Gerrit; Grubmuller, Helmut
  Journal of Chemical Theory and Computation (2011), 7 (6), 1962-1978CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  PH is an important parameter in condensed-phase systems, because it dets. the protonation state of titratable groups and thus influences the structure, dynamics, and function of mols. in soln. In most force field simulation protocols, however, the protonation state of a system (rather than its pH) is kept fixed and cannot adapt to changes of the local environment. Here, we present a method, implemented within the MD package GROMACS, for const. pH mol. dynamics simulations in explicit solvent that is based on the λ-dynamics approach. In the latter, the dynamics of the titrn. coordinate λ, which interpolates between the protonated and deprotonated states, is driven by generalized forces between the protonated and deprotonated states. The hydration free energy, as a function of pH, is included to facilitate const. pH simulations. The protonation states of titratable groups are allowed to change dynamically during a simulation, thus reproducing av. protonation probabilities at a certain pH. The accuracy of the method is tested against titrn. curves of single amino acids and a dipeptide in explicit solvent.
24. 24
  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, DOI: 10.1021/ct200146j
  24
  Continuous Constant pH Molecular Dynamics in Explicit Solvent with pH-Based Replica Exchange
  Wallace, Jason A.; Shen, Jana K.
  Journal of Chemical Theory and Computation (2011), 7 (8), 2617-2629CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  A computational tool that offers accurate pKa values and atomically detailed knowledge of protonation-coupled conformational dynamics is valuable for elucidating mechanisms of energy transduction processes in biol., such as enzyme catalysis and electron transfer as well as proton and drug transport. Toward this goal the authors present a new technique of embedding continuous const. pH mol. dynamics within an explicit-solvent representation. In this technique the authors make use of the efficiency of the generalized-Born (GB) implicit-solvent model for estg. the free energy of protein solvation while propagating conformational dynamics using the more accurate explicit-solvent model. Also, the authors employ a pH-based replica exchange scheme to significantly enhance both protonation and conformational state sampling. Benchmark data of five proteins including HP36, NTL9, BBL, HEWL, and SNase yield an av. abs. deviation of 0.53 and a root mean squared deviation of 0.74 from exptl. data. This level of accuracy is obtained with 1 ns simulations per replica. Detailed anal. reveals that explicit-solvent sampling provides increased accuracy relative to the previous GB-based method by preserving the native structure, providing a more realistic description of conformational flexibility of the hydrophobic cluster, and correctly modeling solvent mediated ion-pair interactions. Thus, the authors anticipate that the new technique will emerge as a practical tool to capture ionization equil. while enabling an intimate view of ionization coupled conformational dynamics that is difficult to delineate with exptl. techniques alone.
25. 25
  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, DOI: 10.1063/1.4766352
  25
  Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH
  Wallace, Jason A.; Shen, Jana K.
  Journal of Chemical Physics (2012), 137 (18), 184105/1-184105/9CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  Recent development of const. pH mol. dynamics (CpHMD) methods has offered promise for adding pH-stat in mol. dynamics simulations. However, until now the working pH mol. dynamics (pHMD) implementations are dependent in part or whole on implicit-solvent models. Here we show that proper treatment of long-range electrostatics and maintaining charge neutrality of the system are crit. for extending the continuous pHMD framework to the all-atom representation. The former is achieved here by adding forces to titrn. coordinates due to long-range electrostatics based on the generalized reaction field method, while the latter is made possible by a charge-leveling technique that couples proton titrn. with simultaneous ionization or neutralization of a co-ion in soln. We test the new method using the pH-replica-exchange CpHMD simulations of a series of aliph. dicarboxylic acids with varying carbon chain length. The av. abs. deviation from the exptl. pKa values is merely 0.18 units. The results show that accounting for the forces due to extended electrostatics removes the large random noise in propagating titrn. coordinates, while maintaining charge neutrality of the system improves the accuracy in the calcd. electrostatic interaction between ionizable sites. Thus, we believe that the way is paved for realizing pH-controlled all-atom mol. dynamics in the near future. (c) 2012 American Institute of Physics.
26. 26
  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, DOI: 10.1021/ct2006314
  26
  Constant pH Molecular Dynamics Simulations of Nucleic Acids in Explicit Solvent
  Goh, Garrett B.; Knight, Jennifer L.; Brooks, Charles L.
  Journal of Chemical Theory and Computation (2012), 8 (1), 36-46CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  The nucleosides of adenine and cytosine have pKa values of 3.50 and 4.08, resp., and are assumed to be unprotonated under physiol. conditions. However, evidence from recent NMR and X-ray crystallog. studies has revealed the prevalence of protonated adenine and cytosine in RNA macromols. Such nucleotides with elevated pKa values may play a role in stabilizing RNA structure and participate in the mechanism of ribozyme catalysis. With the work presented here, we establish the framework and demonstrate the first const. pH MD simulations (CPHMD) for nucleic acids in explicit solvent, in which the protonation state is coupled to the dynamical evolution of the RNA system via λ-dynamics. We adopt the new functional form λNexp for λ that was recently developed for multisite λ-dynamics (MSλD) and demonstrate good sampling characteristics in which rapid and frequent transitions between the protonated and unprotonated states at pH = pKa are achieved. Our calcd. pKa values of simple nucleotides are in a good agreement with exptl. measured values, with a mean abs. error of 0.24 pKa units. This work demonstrates that CPHMD can be used as a powerful tool to investigate pH-dependent biol. properties of RNA macromols.
27. 27
  Chen, W.; Wallace, J. A.; Yue, Z.; Shen, J. K. Introducing Titratable Water to All-Atom Molecular Dynamics at Constant pH. Biophys. J. 2013, 105, L15– L17, DOI: 10.1016/j.bpj.2013.06.036
  27
  Introducing Titratable Water to All-Atom Molecular Dynamics at Constant pH
  Chen, Wei; Wallace, Jason A.; Yue, Zhi; Shen, Jana K.
  Biophysical Journal (2013), 105 (4), L15-L17CODEN: BIOJAU; ISSN:0006-3495. (Cell Press)
  Recent development of titratable coions has paved the way for realizing all-atom mol. dynamics at const. pH. To further improve phys. realism, here we describe a technique in which proton titrn. of the solute is directly coupled to the interconversion between water and hydroxide or hydronium. We test the new method in replica-exchange continuous const. pH mol. dynamics simulations of three proteins, HP36, BBL, and HEWL. The calcd. pKa values based on 10-ns sampling per replica have the av. abs. and root-mean-square errors of 0.7 and 0.9 pH units, resp. Introducing titratable water in mol. dynamics offers a means to model proton exchange between solute and solvent, thus opening a door to gaining new insights into the intricate details of biol. phenomena involving proton translocation.
28. 28
  Lee, J.; Miller, B. T.; Damjanovic, A.; Brooks, B. R. Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange. J. Chem. Theory Comput. 2014, 10, 2738– 2750, DOI: 10.1021/ct500175m
  28
  Constant pH Molecular Dynamics in Explicit Solvent with Enveloping Distribution Sampling and Hamiltonian Exchange
  Lee, Juyong; Miller, Benjamin T.; Damjanovic, Ana; Brooks, Bernard R.
  Journal of Chemical Theory and Computation (2014), 10 (7), 2738-2750CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  We present a new computational approach for const. pH simulations in explicit solvent based on the combination of the enveloping distribution sampling (EDS) and Hamiltonian replica exchange (HREX) methods. Unlike const. pH methods based on variable and continuous charge models, our method is based on discrete protonation states. EDS generates a hybrid Hamiltonian of different protonation states. A smoothness parameter s is used to control the heights of energy barriers of the hybrid-state energy landscape. A small s value facilitates state transitions by lowering energy barriers. Replica exchange between EDS potentials with different s values allows us to readily obtain a thermodynamically accurate ensemble of multiple protonation states with frequent state transitions. The anal. is performed with an ensemble obtained from an EDS Hamiltonian without smoothing, s = ∞, which strictly follows the min. energy surface of the end states. The accuracy and efficiency of this method is tested on aspartic acid, lysine, and glutamic acid, which have two protonation states, a histidine with three states, a four-residue peptide with four states, and snake cardiotoxin with eight states. The pKa values estd. with the EDS-HREX method agree well with the exptl. pKa values. The mean abs. errors of small benchmark systems range from 0.03 to 0.17 pKa units, and those of three titratable groups of snake cardiotoxin range from 0.2 to 1.6 pKa units. This study demonstrates that EDS-HREX is a potent theor. framework, which gives the correct description of multiple protonation states and good calcd. pKa values.
29. 29
  Dobrev, P.; Donnini, S.; Groenhof, G.; Grubmüller, H. Accurate Three States Model for Amino Acids with Two Chemically Coupled Titrating Sites in Explicit Solvent Atomistic Constant pH Simulations and p K a Calculations. J. Chem. Theory Comput. 2017, 13, 147– 160, DOI: 10.1021/acs.jctc.6b00807
  29
  Accurate Three States Model for Amino Acids with Two Chemically Coupled Titrating Sites in Explicit Solvent Atomistic Constant pH Simulations and pKa Calculations
  Dobrev, Plamen; Donnini, Serena; Groenhof, Gerrit; Grubmueller, Helmut
  Journal of Chemical Theory and Computation (2017), 13 (1), 147-160CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  Correct protonation of titratable groups in biomols. is crucial for their accurate description by mol. dynamics simulations. In the context of const. pH simulations, an addnl. protonation degree of freedom is introduced for each titratable site, allowing the protonation state to change dynamically with changing structure or electrostatics. Here, we implement a second reaction coordinate which switches between two tautomeric states of an amino acid with chem. coupled titratable sites, such as aspartate (Asp), glutamate (Glu), and histidine (His). To this aim, we test a scheme involving three protonation states. To facilitate charge neutrality as required for periodic boundary conditions and Particle Mesh Ewald (PME) electrostatics, titrn. of each resp. amino acid is coupled to a "water" mol. that is charged in opposite direction. Addnl., a force field modification for Amber99sb is introduced and tested for the description of carboxyl group protonation. Our three states model is tested by titrn. simulations of Asp, Glu, and His, yielding a good agreement, reproducing the correct geometry of the groups in their different protonation forms. We further show that the ion concn. change due to the neutralizing "water" mols. does not significantly affect the protonation free energies of the titratable groups, suggesting that the three states model provides a good description of biomol. dynamics at const. pH.
30. 30
  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, DOI: 10.1021/acs.jctc.6b00552
  30
  All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water
  Huang, Yandong; Chen, Wei; Wallace, Jason A.; Shen, Jana
  Journal of Chemical Theory and Computation (2016), 12 (11), 5411-5421CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  Development of a pH stat to properly control soln. pH in biomol. simulations has been a long-standing goal in the community. Towards this goal recent years have witnessed the emergence of the so-called const. pH mol. dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or electrostatic truncation schemes. Here the authors report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous const. pH mol. dynamics (CpHMD) method, enabling CpHMD to be performed with a std. MD engine at a fractional added computational cost. The authors demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP38, BBL, HEWL and SNase. With the sampling time of 10 ns per replica, most pKa's are converged, yielding the av. abs. and root-mean-square deviations of 0.61 and 0.77, resp., from expt. Linear regression of the calcd. vs. exptl. pKa shifts gives a correlation coeff. of 0.79, a slope of 1 and an intercept near 0. Anal. reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromol. systems in various environments, enabling at. insights into pH-dependent phenomena involving not only sol. proteins but also transmembrane proteins, nucleic acids, surfactants and polysaccharides.
31. 31
  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, DOI: 10.1021/acs.jcim.8b00227
  31
  Generalized Born Based Continuous Constant pH Molecular Dynamics in Amber: Implementation, Benchmarking and Analysis
  Huang, Yandong; Harris, Robert C.; Shen, Jana
  Journal of Chemical Information and Modeling (2018), 58 (7), 1372-1383CODEN: JCISD8; ISSN:1549-9596. (American Chemical Society)
  Soln. pH plays an important role in structure and dynamics of biomol. systems; however, pH effects cannot be accurately accounted for in conventional mol. dynamics simulations based on fixed protonation states. Continuous const. pH mol. dynamics (CpHMD) based on the λ-dynamics framework calcs. protonation states on the fly during dynamical simulation at a specified pH condition. Here the authors report the CPU-based implementation of the CpHMD method based on the GBNeck2 generalized Born (GB) implicit-solvent model in the pmemd engine of the Amber mol. dynamics package. The performance of the method was tested using pH replica-exchange titrn. simulations of Asp, Glu and His side chains in 4 miniproteins and 7 enzymes with exptl. known pKa's, some of which are significantly shifted from the model values. The added computational cost due to CpHMD titrn. ranges from 11 to 33% for the data set and scales roughly linearly as the ratio between the titrable sites and no. of solute atoms. Comparison of the exptl. and calcd. pKa's using 2 ns per replica sampling yielded a mean unsigned error of 0.70, a root-mean-squared error of 0.91, and a linear correlation coeff. of 0.79. Though this level of accuracy is similar to the GBSW-based CpHMD in CHARMM, in contrast to the latter, the current implementation was able to reproduce the exptl. orders of the pKa's of the coupled carboxylic dyads. The authors quantified the sampling errors, which revealed that prolonged simulation is needed to converge pKa's of several titratable groups involved in salt-bridge-like interactions or deeply buried in the protein interior. The authors' benchmark data demonstrate that GBNeck2-CpHMD is an attractive tool for protein pKa predictions.
32. 32
  Harris, R. C.; Shen, J. GPU-accelerated implementation of continuous constant pH molecular dynamics in amber: pKa predictions with single-pH simulations. J. Chem. Inf. Model. 2019, 59, 4821– 4832, DOI: 10.1021/acs.jcim.9b00754
  32
  GPU-Accelerated Implementation of Continuous Constant pH Molecular Dynamics in Amber: pKa Predictions with Single-pH Simulations
  Harris, Robert C.; Shen, Jana
  Journal of Chemical Information and Modeling (2019), 59 (11), 4821-4832CODEN: JCISD8; ISSN:1549-9596. (American Chemical Society)
  The authors present a GPU implementation of the continuous const. pH mol. dynamics (CpHMD) based on the most recent generalized Born implicit-solvent model in the pmemd engine of the Amber mol. dynamics package. To test the accuracy of the tool for rapid pKa predictions, a series of 2-ns single-pH simulations were performed for over 120 titratable residues in 10 benchmark proteins that were previously used to test the various continuous CpHMD methods. The calcd. pKa's showed a root-mean-square deviation of 0.80 and correlation coeff. of 0.83 with respect to expt. 90% of the pKa's were converged with estd. errors <0.1 pH units. Surprisingly, this level of accuracy is similar to the authors' previous replica-exchange simulations with 2 ns per replica and an exchange attempt frequency of 2 ps-1 (Huang, Harris and Shen, J Chem Info Model, 2018). For the linked titrn. sites in two enzymes, although residue-specific protonation state sampling in the single-pH simulations was not converged within 2 ns, the protonation fraction of the linked residues appeared to be largely converged, and the exptl. macroscopic {\pka} values were reproduced to within 1 pH unit. Comparison with replica-exchange simulations with different exchange attempt frequencies showed that the splitting between the two macroscopic pKa's is underestimated with frequent exchange attempts such as 2 ps$̂{-1}$, while single-pH simulations overestimate the splitting. The same trend is seen for the single-pH vs. replica-exchange simulations of a hydrogen-bonded aspartyl dyad in a much larger protein. A 2-ns single-pH simulation of a 400-residue protein takes about one hour on a single NVIDIA GeForce RTX 2080 graphics card, which is over 1000 times faster than a CpHMD run on a single CPU core of a high-performance computing cluster node. Thus, the authors envision that GPU-accelerated continuous CpHMD may be used in routine pKa predictions for a variety of applications, from assisting MD simulations with protonation state assignment to offering pH-dependent corrections of binding free energies and identifying reactive hot spots for covalent drug design.
33. 33
  Harris, R. C.; Liu, R.; Shen, J. Predicting reactive cysteines with implicit-solvent-based continuous constant pH molecular dynamics in amber. J. Chem. Theory Comput. 2020, 16, 3689– 3698, DOI: 10.1021/acs.jctc.0c00258
  33
  Predicting Reactive Cysteines with Implicit-Solvent-Based Continuous Constant pH Molecular Dynamics in Amber
  Harris, Robert C.; Liu, Ruibin; Shen, Jana
  Journal of Chemical Theory and Computation (2020), 16 (6), 3689-3698CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  Cysteines existing in the deprotonated thiolate form or having a tendency to become deprotonated are important players in enzymic and cellular redox functions and frequently exploited in covalent drug design; however, most computational studies assume cysteines as protonated. Thus, developing an efficient tool that can make accurate and reliable predictions of cysteine protonation states is timely needed. The authors recently implemented a generalized Born (GB) based continuous const. pH mol. dynamics (CpHMD) method in Amber for protein pKa calcns. on CPUs and GPUs. Here the authors benchmark the performance of GB-CpHMD for predictions of cysteine pKa's and reactivities using a data set of 24 proteins with both down- and upshifted cysteine pKa's. The authors found that 10 ns single-pH or 4 ns replica-exchange CpHMD titrns. gave root-mean-square errors of 1.2-1.3 and correlation coeffs. of 0.8-0.9 with respect to expt. The accuracy of predicting thiolates or reactive cysteines at physiol. pH with single-pH titrns. is 86 or 81% with a precision of 100 or 90%, resp. This performance well surpasses the traditional structure-based methods, particularly a widely used empirical pKa tool that gives an accuracy less than 50%. The authors discuss simulation convergence, dependence on starting structures, common determinants of the pKa downshifts and upshifts, and the origin of the discrepancies from the structure-based calcns. The work suggests that CpHMD titrns. can be performed on a desktop computer equipped with a single GPU card to predict cysteine protonation states for a variety of applications, from understanding biol. functions to covalent drug design.
34. 34
  Grünewald, F.; Souza, P. C.; Abdizadeh, H.; Barnoud, J.; de Vries, A. H.; Marrink, S. J. Titratable Martini model for constant pH simulations. J. Chem. Phys. 2020, 153, 024118, DOI: 10.1063/5.0014258
  34
  Titratable Martini model for constant pH simulations
  Gruenewald, Fabian; Souza, Paulo C. T.; Abdizadeh, Haleh; Barnoud, Jonathan; de Vries, Alex H.; Marrink, Siewert J.
  Journal of Chemical Physics (2020), 153 (2), 024118CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  The authors deliver a proof of concept for a fast method that introduces pH effects into classical coarse-grained (CG) mol. dynamics simulations. The authors' approach is based upon the latest version of the popular Martini CG model to which explicit proton mimicking particles were added. The authors verify the authors' approach against exptl. data involving several different mols. and different environmental conditions. In particular, the authors compute titrn. curves, pH dependent free energies of transfer, and lipid bilayer membrane affinities as a function of pH. Using oleic acid as an example compd., the authors further illustrate that the authors' method can be used to study passive translocation in lipid bilayers via protonation. Finally, the authors' model reproduces qual. the expansion of the macromol. dendrimer poly(propylene imine) as well as the assocd. pKa shift of its different generations. This example demonstrates that the authors' model is able to pick up collective interactions between titratable sites in large mols. comprising many titratable functional groups. (c) 2020 American Institute of Physics.
35. 35
  Mongan, J.; Case, D. A. Biomolecular simulations at constant pH. Curr. Opin. Struct. Biol. 2005, 15, 157– 163, DOI: 10.1016/j.sbi.2005.02.002
  35
  Biomolecular simulations at constant pH
  Mongan, John; Case, David A.
  Current Opinion in Structural Biology (2005), 15 (2), 157-163CODEN: COSBEF; ISSN:0959-440X. (Elsevier Ltd.)
  A review. Like temp. and pressure, the soln. pH is an important thermodn. variable that is commonly varied in expts. and is used by cells to influence biochem. function. It is now becoming feasible to carry out practical mol. dynamics simulations that mimic the thermodn. of such expts., by allowing proton transfer between the system of interest and a hypothetical bath of protons at a given pH. These are demanding calcns., because the energetics of charge changes upon protonation or deprotonation must be accurately modeled, and because such simulations must sample both mol. configurations and the large no. of protonation states that are possible for a mol. with many titrating sites.
36. 36
  Damjanovic, A.; Miller, B. T.; Okur, A.; Brooks, B. R. Reservoir pH replica exchange. J. Chem. Phys. 2018, 149, 072321, DOI: 10.1063/1.5027413
  36
  Reservoir pH replica exchange
  Damjanovic, Ana; Miller, Benjamin T.; Okur, Asim; Brooks, Bernard R.
  Journal of Chemical Physics (2018), 149 (7), 072321/1-072321/12CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  We present the reservoir pH replica exchange (R-pH-REM) method for const. pH simulations. The R-pH-REM method consists of a two-step procedure; the first step involves generation of one or more reservoirs of conformations. Each reservoir is obtained from a std. or enhanced mol. dynamics simulation with a constrained (fixed) protonation state. In the second step, fixed charge constraints are relaxed, as the structures from one or more reservoirs are periodically injected into a const. pH or a pH-replica exchange (pH-REM) simulation. The benefit of this two-step process is that the computationally intensive part of conformational search can be decoupled from const. pH simulations, and various techniques for enhanced conformational sampling can be applied without the need to integrate such techniques into the pH-REM framework. Simulations on blocked Lys, KK, and KAAE peptides were used to demonstrate an agreement between pH-REM and R-pH-REM simulations. While the reservoir simulations are not needed for these small test systems, the real need arises in cases when ionizable mols. can sample two or more conformations sepd. by a large energy barrier, such that adequate sampling is not achieved on a time scale of std. const. pH simulations. Such problems might be encountered in protein systems that exploit conformational transitions for function. A hypothetical case is studied, a small mol. with a large torsional barrier; while results of pH-REM simulations depend on the starting structure, R-pH-REM calcns. on this model system are in excellent agreement with a theor. model. (c) 2018 American Institute of Physics.
37. 37
  Chen, Y.; Roux, B. Constant-pH Hybrid Nonequilibrium Molecular Dynamics - Monte Carlo Simulation Method. J. Chem. Theory Comput. 2015, 11, 3919– 3931, DOI: 10.1021/acs.jctc.5b00261
  37
  Constant-pH Hybrid Nonequilibrium Molecular Dynamics-Monte Carlo Simulation Method
  Chen, Yunjie; Roux, Benoit
  Journal of Chemical Theory and Computation (2015), 11 (8), 3919-3931CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  A computational method is developed to carry out explicit solvent simulations of complex mol. systems under conditions of const. pH. In const.-pH simulations, preidentified ionizable sites are allowed to spontaneously protonate and deprotonate as a function of time in response to the environment and the imposed pH. The method consists of carrying out short nonequil. mol. dynamics (neMD) switching trajectories to generate phys. plausible configurations with changed protonation states that are subsequently accepted or rejected according to a Metropolis Monte Carlo (MC) criterion. To ensure microscopic detailed balance arising from such nonequil. switches, the at. momenta are altered according to the sym. two-ends momentum reversal prescription. To achieve higher efficiency, the original neMD-MC scheme is sepd. into two steps, reducing the need for generating a large no. of unproductive and costly nonequil. trajectories. In the first step, the protonation state of a site is randomly attributed via a Metropolis MC process on the basis of an intrinsic pKa; an attempted nonequil. switch is generated only if this change in protonation state is accepted. This hybrid two-step inherent pKa neMD-MC simulation method is tested with single amino acids in soln. (Asp, Glu, and His) and then applied to turkey ovomucoid third domain and hen egg-white lysozyme. Because of the simple linear increase in the computational cost relative to the no. of titratable sites, the present method is naturally able to treat extremely large systems.
38. 38
  Kong, X.; Brooks, C. L., III λ-dynamics: A new approach to free energy calculations. J. Chem. Phys. 1996, 105, 2414– 2423, DOI: 10.1063/1.472109
  There is no corresponding record for this reference.
39. 39
  Simonson, T.; Carlsson, J.; Case, D. A. Proton binding to proteins: p K a calculations with explicit and implicit solvent models. J. Am. Chem. Soc. 2004, 126, 4167– 4180, DOI: 10.1021/ja039788m
  39
  Proton Binding to Proteins: pKa Calculations with Explicit and Implicit Solvent Models
  Simonson, Thomas; Carlsson, Jens; Case, David A.
  Journal of the American Chemical Society (2004), 126 (13), 4167-4180CODEN: JACSAT; ISSN:0002-7863. (American Chemical Society)
  Ionizable residues play important roles in protein structure and activity, and proton binding is a valuable reporter of electrostatic interactions in these systems. The authors use mol. dynamics free energy simulations (MDFE) to compute proton pKa shifts, relative to a model compd. in soln., for three aspartate side chains in two proteins. Simulations with explicit solvent and with an implicit, dielec. continuum solvent are reported. The implicit solvent simulations use the generalized Born (GB) model, which provides an approx., anal. soln. to Poisson's equation. With explicit solvent, the direction of the pKa shifts is correct in all three cases with one force field (AMBER) and in two out of three cases with another (CHARMM). For two aspartates, the dielec. response to ionization is found to be linear, even though the sep. protein and solvent responses can be nonlinear. For thioredoxin Asp-26, nonlinearity arises from the presence of two substates that correspond to the two possible orientations of the protonated carboxylate. For this side chain, which is partly buried and has a large pKa upshift, very long simulations are needed to correctly sample several slow degrees of freedom that reorganize in response to the ionization. Thus, nearby Lys-57 rotates to form a salt bridge and becomes buried, while three waters intercalate along the opposite edge of Asp-26. Such strong and anisotropic reorganization is very difficult to predict with Poisson-Boltzmann methods that only consider electrostatic interactions and employ a single protein structure. In contrast, MDFE with a GB dielec. continuum solvent, used for the first time for pKa calcns., can describe protein reorganization accurately and gives encouraging agreement with expt. and with the explicit solvent simulations.
40. 40
  Chen, W.; Shen, J. K. Effects of System Net Charge and Electrostic Truncation on All-Atom Constant pH Molecular Dynamics. J. Comput. Chem. 2014, 35, 1986– 1996, DOI: 10.1002/jcc.23713
  40
  Effects of system net charge and electrostatic truncation on all-atom constant pH molecular dynamics
  Chen, Wei; Shen, Jana K.
  Journal of Computational Chemistry (2014), 35 (27), 1986-1996CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
  Const. pH mol. dynamics offers a means to rigorously study the effects of soln. pH on dynamical processes. Here, we address two crit. questions arising from the most recent developments of the all-atom continuous const. pH mol. dynamics (CpHMD) method: (1) What is the effect of spatial electrostatic truncation on the sampling of protonation states. (2) Is the enforcement of elec. neutrality necessary for const. pH simulations. We first examd. how the generalized reaction field and force-shifting schemes modify the electrostatic forces on the titrn. coordinates. Free energy simulations of model compds. were then carried out to delineate the errors in the deprotonation free energy and salt-bridge stability due to electrostatic truncation and system net charge. Finally, CpHMD titrn. of a mini-protein HP36 was used to understand the manifestation of the two types of errors in the calcd. pKa values. The major finding is that enforcing charge neutrality under all pH conditions and at all time via co-titrating ions significantly improves the accuracy of protonation-state sampling. We suggest that such finding is also relevant for simulations with particle mesh Ewald, considering the known artifacts due to charge-compensating background plasma.
41. 41
  Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; De Vries, A. H. The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812– 7824, DOI: 10.1021/jp071097f
  41
  The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations
  Marrink, Siewert J.; Risselada, H. Jelger; Yefimov, Serge; Tieleman, D. Peter; De Vries, Alex H.
  Journal of Physical Chemistry B (2007), 111 (27), 7812-7824CODEN: JPCBFK; ISSN:1520-6106. (American Chemical Society)
  We present an improved and extended version of our coarse grained lipid model. The new version, coined the MARTINI force field, is parametrized in a systematic way, based on the reprodn. of partitioning free energies between polar and apolar phases of a large no. of chem. compds. To reproduce the free energies of these chem. building blocks, the no. of possible interaction levels of the coarse-grained sites has increased compared to those of the previous model. Application of the new model to lipid bilayers shows an improved behavior in terms of the stress profile across the bilayer and the tendency to form pores. An extension of the force field now also allows the simulation of planar (ring) compds., including sterols. Application to a bilayer/cholesterol system at various concns. shows the typical cholesterol condensation effect similar to that obsd. in all atom representations.
42. 42
  Donnini, S.; Ullmann, R. T.; Groenhof, G.; Grubmüller, H. Charge-neutral constant pH molecular dynamics simulations using a parsimonious proton buffer. J. Chem. Theory Comput. 2016, 12, 1040– 1051, DOI: 10.1021/acs.jctc.5b01160
  42
  Charge-Neutral Constant pH Molecular Dynamics Simulations Using a Parsimonious Proton Buffer
  Donnini, Serena; Ullmann, R. Thomas; Groenhof, Gerrit; Grubmuller, Helmut
  Journal of Chemical Theory and Computation (2016), 12 (3), 1040-1051CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  In const. pH mol. dynamics simulations, the protonation states of titratable sites can respond to changes of the pH and of their electrostatic environment. Consequently, the no. of protons bound to the biomol., and therefore the overall charge of the system, fluctuates during the simulation. To avoid artifacts assocd. with a non-neutral simulation system, we introduce an approach to maintain neutrality of the simulation box in const. pH mol. dynamics simulations, while maintaining an accurate description of all protonation fluctuations. Specifically, we introduce a proton buffer that, like a buffer in expt., can exchange protons with the biomol. enabling its charge to fluctuate. To keep the total charge of the system const., the uptake and release of protons by the buffer are coupled to the titrn. of the biomol. with a constraint. We find that, because the fluctuation of the total charge (no. of protons) of a typical biomol. is much smaller than the no. of titratable sites of the biomol., the no. of buffer sites required to maintain overall charge neutrality without compromising the charge fluctuations of the biomol., is typically much smaller than the no. of titratable sites, implying markedly enhanced simulation and sampling efficiency.
43. 43
  Dobrev, P.; Vemulapalli, S. P. B.; Nath, N.; Griesinger, C.; Grubmüller, H. Probing the accuracy of explicit solvent constant pH molecular dynamics simulations for peptides. J. Chem. Theory Comput. 2020, 16, 2561– 2569, DOI: 10.1021/acs.jctc.9b01232
  43
  Probing the Accuracy of Explicit Solvent Constant pH Molecular Dynamics Simulations for Peptides
  Dobrev, Plamen; Vemulapalli, Sahithya Phani Babu; Nath, Nilamoni; Griesinger, Christian; Grubmueller, Helmut
  Journal of Chemical Theory and Computation (2020), 16 (4), 2561-2569CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  Protonation states of titratable amino acids play a key role in many biomol. processes. Knowledge of protonatable residue charges at a given pH is essential for a correct understanding of protein catalysis, inter- and intramol. interactions, substrate binding, and protein dynamics for instance. However, acquiring exptl. values for individual amino acid protonation states of complex systems is not straightforward; therefore, several in silico approaches have been developed to tackle this issue. In this work, the authors assess the accuracy of the previously developed const. pH MD approach by comparing the theor. obtained pKa values for titratable residues with exptl. values from an equiv. NMR study. The authors selected a set of four pentapeptides, of adequately small size to ensure comprehensive sampling, but concurrently, due to their charge compn., posing a challenge for protonation state calcn. The comparison of the pKa values shows good agreement of the exptl. and the theor. approach with a largest difference of 0.25 pKa units. Further, the corresponding titrn. curves are in fair agreement, although the shift of the Hill coeff. from a value of 1 was not always reproduced in simulations. The phase space overlap in Cartesian space between trajectories generated in const. pH and std. MD simulations is fair and suggests that the const. pH MD approach reasonably well preserves the dynamics of the system, allowing dynamic protonation MD simulations without introducing structural artifacts.
44. 44
  Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089– 10092, DOI: 10.1063/1.464397
  44
  Particle mesh Ewald: an N·log(N) method for Ewald sums in large systems
  Darden, Tom; York, Darrin; Pedersen, Lee
  Journal of Chemical Physics (1993), 98 (12), 10089-92CODEN: JCPSA6; ISSN:0021-9606.
  An N·log(N) method for evaluating electrostatic energies and forces of large periodic systems is presented. The method is based on interpolation of the reciprocal space Ewald sums and evaluation of the resulting convolution using fast Fourier transforms. Timings and accuracies are presented for three large cryst. ionic systems.
45. 45
  Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577– 8593, DOI: 10.1063/1.470117
  45
  A smooth particle mesh Ewald method
  Essmann, Ulrich; Perera, Lalith; Berkowitz, Max L.; Darden, Tom; Lee, Hsing; Pedersen, Lee G.
  Journal of Chemical Physics (1995), 103 (19), 8577-93CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  The previously developed particle mesh Ewald method is reformulated in terms of efficient B-spline interpolation of the structure factors. This reformulation allows a natural extension of the method to potentials of the form 1/rp with p ≥ 1. Furthermore, efficient calcn. of the virial tensor follows. Use of B-splines in the place of Lagrange interpolation leads to analytic gradients as well as a significant improvement in the accuracy. The authors demonstrate that arbitrary accuracy can be achieved, independent of system size N, at a cost that scales as N log(N). For biomol. systems with many thousands of atoms and this method permits the use of Ewald summation at a computational cost comparable to that of a simple truncation method of 10 Å or less.
46. 46
  Knight, J. L.; Brooks, C. L., III Multisite λ dynamics for simulated structure-activity relationship studies. J. Chem. Theory Comput. 2011, 7, 2728– 2739, DOI: 10.1021/ct200444f
  46
  Multisite λ Dynamics for Simulated Structure-Activity Relationship Studies
  Knight, Jennifer L.; Brooks, Charles L., III
  Journal of Chemical Theory and Computation (2011), 7 (9), 2728-2739CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  Multisite λ dynamics (MSλD) is a new free energy simulation method that is based on λ dynamics. It has been developed to enable multiple substituents at multiple sites on a common ligand core to be modeled simultaneously and their free energies assessed. The efficacy of MSλD for estg. relative hydration free energies and relative binding affinties is demonstrated using three test systems. Model compds. representing multiple identical benzene, dihydroxybenzene, and dimethoxybenzene mols. show that total combined MSλD trajectory lengths of ∼1.5 ns are sufficient to reliably achieve relative hydration free energy ests. within 0.2 kcal/mol and are less sensitive to the no. of trajectories that are used to generate these ests. for hybrid ligands that contain up to 10 substituents modeled at a single site or five substituents modeled at each of two sites. Relative hydration free energies among six benzene derivs. calcd. from MSλD simulations are in very good agreement with those from alchem. free energy simulations (with av. unsigned differences of 0.23 kcal/mol and R2 = 0.991) and the expt. (with av. unsigned errors of 1.8 kcal/mol and R2 = 0.959). Ests. of the relative binding affinities among 14 inhibitors of HIV-1 reverse transcriptase obtained from MSλD simulations are in reasonable agreement with those from traditional free energy simulations and the expt. (av. unsigned errors of 0.9 kcal/mol and R2 = 0.402). For the same level of accuracy and precision, MSλD simulations are achieved ∼20-50 times faster than traditional free energy simulations and thus with reliable force field parameters can be used effectively to screen tens to hundreds of compds. in structure-based drug design applications.
47. 47
  Goh, G. B.; Hulbert, B. S.; Zhou, H.; Brooks, C. L., III Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism. Proteins: Struct., Funct., Bioinf. 2014, 82, 1319– 1331, DOI: 10.1002/prot.24499
  47
  Constant pH molecular dynamics of proteins in explicit solvent with proton tautomerism
  Goh, Garrett B.; Hulbert, Benjamin S.; Zhou, Huiqing; Brooks, Charles L., III
  Proteins: Structure, Function, and Bioinformatics (2014), 82 (7), 1319-1331CODEN: PSFBAF ISSN:. (Wiley-Blackwell)
  PH is a ubiquitous regulator of biol. activity, including protein-folding, protein-protein interactions, and enzymic activity. Existing const. pH mol. dynamics (CPHMD) models that were developed to address questions related to the pH-dependent properties of proteins are largely based on implicit solvent models. However, implicit solvent models are known to underestimate the desolvation energy of buried charged residues, increasing the error assocd. with predictions that involve internal ionizable residue that are important in processes like hydrogen transport and electron transfer. Furthermore, discrete water and ions cannot be modeled in implicit solvent, which are important in systems like membrane proteins and ion channels. We report on an explicit solvent const. pH mol. dynamics framework based on multi-site λ-dynamics (CPHMDMSλD). In the CPHMDMSλD framework, we performed seamless alchem. transitions between protonation and tautomeric states using multi-site λ-dynamics, and designed novel biasing potentials to ensure that the phys. end-states are predominantly sampled. We show that explicit solvent CPHMDMSλD simulations model realistic pH-dependent properties of proteins such as the Hen-Egg White Lysozyme (HEWL), binding domain of 2-oxoglutarate dehydrogenase (BBL) and N-terminal domain of ribosomal protein L9 (NTL9), and the pKa predictions are in excellent agreement with exptl. values, with a RMSE ranging from 0.72 to 0.84 pKa units. With the recent development of the explicit solvent CPHMDMSλD framework for nucleic acids, accurate modeling of pH-dependent properties of both major class of biomols.-proteins and nucleic acids is now possible.
48. 48
  Knight, J. L.; Brooks, C. L., III Applying efficient implicit nongeometric constraints in alchemical free energy simulations. Journal of computational chemistry 2011, 32, 3423– 3432, DOI: 10.1002/jcc.21921
  48
  Applying efficient implicit nongeometric constraints in alchemical free energy simulations
  Knight, Jennifer L.; Brooks, Charles L.
  Journal of Computational Chemistry (2011), 32 (16), 3423-3432CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
  Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in mol. dynamics simulations. Advanced variations of alchem. free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are assocd. with the progress variables in alchem. transformations are treated as dynamic variables and are constrained such that: 0 ≤ λi ≤ 1 and .sum.i λi = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by a functional form of λi that oscillates between λi = 0 and λi = 1 and has relatively steep transitions between these endpoints. Copyright for JCC Journal: © 2011 Wiley Periodicals, Inc.; J. Comput. Chem., 2011.
49. 49
  Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015, 1, 19– 25, DOI: 10.1016/j.softx.2015.06.001
  There is no corresponding record for this reference.
50. 50
  Singhal, A. K.; Chien, K. Y.; Wu, W. G.; Rule, G. S. Solution structure of cardiotoxin V from Naja naja atra. Biochemistry 1993, 32, 8036– 8044, DOI: 10.1021/bi00082a026
  50
  Solution structure of cardiotoxin V from Naja naja atra
  Singhal, Arun K.; Chien, Kun Yi; Wu, Wen Guey; Rule, Gordon S.
  Biochemistry (1993), 32 (31), 8036-44CODEN: BICHAW; ISSN:0006-2960.
  To det. whether the unique activity of this cardiotoxin (CTX) is attributable to its tertiary structure, the soln. structure of CTX V was detd. by NMR methods. On the basis of these studies, this cardiotoxin has the same general topol. as other members of the family, and thus its unusual properties do not arise from any gross structural differences that are detectable by soln. NMR methods. Mol. dynamics calcns. indicate that residues 36-50 show concerted fluctuations. On the basis of sequence similarity, the authors postulate that residues 30-34 are important in detg. the specificity of cardiotoxins for fusion vs. lysis of vesicles.
51. 51
  Ramanadham, M.; Sieker, L.; Jensen, L. Refinement of triclinic lysozyme: II. The method of stereochemically restrained least squares. Acta Crystallographica Section B: Structural Science 1990, 46, 63– 69, DOI: 10.1107/S0108768189009195
  51
  Refinement of triclinic lysozyme: II. The method of stereochemically restrained least squares
  Ramanadham M; Sieker L C; Jensen L H
  Acta crystallographica. Section B, Structural science (1990), 46 ( Pt 1) (), 63-9 ISSN:0108-7681.
  Refinement of triclinic lysozyme by restrained least squares against the 2 A resolution X-ray data is described, beginning with the model from cycle 17 of the preceding paper [Hodsdon, Brown, Sieker & Jensen (1990). Acta Cryst. B46, 54-62]. After 20 refinement cycles, R stood at 0.172. Nevertheless, serious errors involving both main-chain and side-chain atoms still remained, requiring numerous model rebuilding sessions interleaved with refinement cycles. After 63 cycles R = 0.124 for the model which includes all protein atoms, 249 water oxygen sites and five nitrate ions. Although the overall B is relatively low, 10.5 A2, B's for atoms in the region of residues 101-103, toward the termini of some of the longer side chains, and in the region of the C terminus of the main chain exceed 20 A2, indicating relatively high atomic mobilities, disorder, or remaining errors in the model.
52. 52
  Sauguet, L.; Poitevin, F.; Murail, S.; Van Renterghem, C.; Moraga-Cid, G.; Malherbe, L.; Thompson, A. W.; Koehl, P.; Corringer, P.-J.; Baaden, M.; Delarue, M. Structural basis for ion permeation mechanism in pentameric ligand-gated ion channels. EMBO journal 2013, 32, 728– 741, DOI: 10.1038/emboj.2013.17
  52
  Structural basis for ion permeation mechanism in pentameric ligand-gated ion channels
  Sauguet, Ludovic; Poitevin, Frederic; Murail, Samuel; Van Renterghem, Catherine; Moraga-Cid, Gustavo; Malherbe, Laurie; Thompson, Andrew W.; Koehl, Patrice; Corringer, Pierre-Jean; Baaden, Marc; Delarue, Marc
  EMBO Journal (2013), 32 (5), 728-741CODEN: EMJODG; ISSN:0261-4189. (Nature Publishing Group)
  To understand the mol. mechanism of ion permeation in pentameric ligand-gated ion channels (pLGIC), the structure of an open form of GLIC, a prokaryotic pLGIC, was solved at 2.4 Å. Anomalous diffraction data were used to place bound anions and cations. This reveals ordered water mols. at the level of two rings of hydroxylated residues (named Ser6' and Thr2') that contribute to the ion selectivity filter. Two water pentagons are obsd., a self-stabilized ice-like water pentagon and a second wider water pentagon, with one sodium ion between them. Single-channel electrophysiol. shows that the side-chain hydroxyl of Ser6' is crucial for ion translocation. Simulations and electrostatics calcns. complemented the description of hydration in the pore and suggest that the water pentagons obsd. in the crystal are important for the ion to cross hydrophobic constriction barriers. Simulations that pull a cation through the pore reveal that residue Ser6' actively contributes to ion translocation by reorienting its side chain when the ion is going through the pore. Generalization of these findings to the pLGIC family is proposed.
53. 53
  Lee, T.-W.; Qasim, M., Jr; Laskowski, M., Jr; James, M. N. Structural insights into the non-additivity effects in the sequence-to-reactivity algorithm for serine peptidases and their inhibitors. Journal of molecular biology 2007, 367, 527– 546, DOI: 10.1016/j.jmb.2007.01.008
  53
  Structural Insights into the Non-additivity Effects in the Sequence-to-Reactivity Algorithm for Serine Peptidases and their Inhibitors
  Lee, Ting-Wai; Qasim, M. A.; Laskowski, Michael; James, Michael N. G.
  Journal of Molecular Biology (2007), 367 (2), 527-546CODEN: JMOBAK; ISSN:0022-2836. (Elsevier Ltd.)
  Sequence-to-reactivity algorithms (SRAs) for proteins have the potential of being broadly applied in mol. design. Recently, Laskowski et al. have reported an additivity-based SRA that accurately predicts most of the std. free energy changes of assocn. for variants of turkey ovomucoid third domain (OMTKY3) with six serine peptidases, one of which is streptogrisin B (commonly known as Streptomyces griseus peptidase B, SGPB). Non-additivity effects for residues 18I and 32I, and for residues 20I and 32I of OMTKY3 occurred when the assocns. with SGPB were predicted using the SRA. To elucidate precisely the mechanics of these non-additivity effects in structural terms, we have detd. the crystal structures of the unbound OMTKY3 (with Gly32I as in the wild-type amino acid sequence) at a resoln. of 1.16 Å, the unbound Ala32I variant of OMTKY3 at a resoln. of 1.23 Å, and the SGPB:OMTKY3-Ala32I complex (equil. assocn. const. Ka = 7.1×109 M-1 at 21(±2) C°, pH 8.3) at a resoln. of 1.70 Å. Extensive comparisons with the crystal structure of the unbound OMTKY3 confirm our understanding of some previously addressed non-additivity effects. Unexpectedly, SGPB and OMTKY3-Ala32I form a 1:2 complex in the crystal. Comparison with the SGPB:OMTKY3 complex shows a conformational change in the SGPB:OMTKY3-Ala32I complex, resulting from a hinged rigid-body rotation of the inhibitor caused by the steric hindrance between the Me group of Ala32IA of the inhibitor and Pro192BE of the peptidase. This perturbs the interactions among residues 18I, 20I, 32I and 36I of the inhibitor, probably resulting in the above non-additivity effects. This conformational change also introduces residue 10I as an addnl. hyper-variable contact residue to the SRA.
54. 54
  Huang, J.; MacKerell, A. D., Jr CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data. Journal of computational chemistry 2013, 34, 2135– 2145, DOI: 10.1002/jcc.23354
  54
  CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data
  Huang, Jing; MacKerell, Alexander D. Jr
  Journal of Computational Chemistry (2013), 34 (25), 2135-2145CODEN: JCCHDD; ISSN:0192-8651. (John Wiley & Sons, Inc.)
  Protein structure and dynamics can be characterized on the atomistic level with both NMR expts. and mol. dynamics (MD) simulations. Here, the authors quantify the ability of the recently presented CHARMM36 (C36) force field (FF) to reproduce various NMR observables using MD simulations. The studied NMR properties include backbone scalar couplings across hydrogen bonds, residual dipolar couplings (RDCs) and relaxation order parameter, as well as scalar couplings, RDCs, and order parameters for side-chain amino- and methyl-contg. groups. The C36 FF leads to better correlation with exptl. data compared to the CHARMM22/CMAP FF and suggest using C36 in protein simulations. Although both CHARMM FFs contains the same nonbond parameters, the authors' results show how the changes in the internal parameters assocd. with the peptide backbone via CMAP and the χ1 and χ2 dihedral parameters leads to improved treatment of the analyzed nonbond interactions. This highlights the importance of proper treatment of the internal covalent components in modeling nonbond interactions with mol. mechanics FFs.
55. 55
  Buslaev, P.; Aho, N.; Jansen, A.; Bauer, P.; Hess, B.; Groenhof, G. Best practices in constant pH MD simulations: accuracy and precision. J. Chem. Theory Comput. 2022, DOI: 10.1021/acs.jctc.2c00517 .
  There is no corresponding record for this reference.
56. 56
  Martini-website, General purpose coarse-grained force field. http://cgmartini.nl/index.php/tools2/proteins-and-bilayers/204-martinize (accessed October 01, 2021).
  There is no corresponding record for this reference.
57. 57
  Kaminski, G. A.; Friesner, R. A.; Tirado-Rives, J.; Jorgensen, W. L. Evaluation and reparametrization of the OPLS-AA force field for proteins via comparison with accurate quantum chemical calculations on peptides. J. Phys. Chem. B 2001, 105, 6474– 6487, DOI: 10.1021/jp003919d
  57
  Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides
  Kaminski, George A.; Friesner, Richard A.; Tirado-Rives, Julian; Jorgensen, William L.
  Journal of Physical Chemistry B (2001), 105 (28), 6474-6487CODEN: JPCBFK; ISSN:1089-5647. (American Chemical Society)
  We present results of improving the OPLS-AA force field for peptides by means of refitting the key Fourier torsional coeffs. The fitting technique combines using accurate ab initio data as the target, choosing an efficient fitting subspace of the whole potential-energy surface, and detg. wts. for each of the fitting points based on magnitudes of the potential-energy gradient. The av. energy RMS deviation from the LMP2/cc-pVTZ(-f)//HF/6-31G** data is reduced by ∼40% from 0.81 to 0.47 kcal/mol as a result of the fitting for the electrostatically uncharged dipeptides. Transferability of the parameters is demonstrated by using the same alanine dipeptide-fitted backbone torsional parameters for all of the other dipeptides (with the appropriate side-chain refitting) and the alanine tetrapeptide. Parameters of nonbonded interactions have also been refitted for the sulfur-contg. dipeptides (cysteine and methionine), and the validity of the new Coulombic charges and the van der Waals σ's and ε's is proved through reproducing gas-phase energies of complex formation heats of vaporization and densities of pure model liqs. Moreover, a novel approach to fitting torsional parameters for electrostatically charged mol. systems has been presented and successfully tested on five dipeptides with charged side chains.
58. 58
  Buck, M.; Bouguet-Bonnet, S.; Pastor, R. W.; MacKerell, A. D., Jr Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme. Biophysical journal 2006, 90, L36– L38, DOI: 10.1529/biophysj.105.078154
  58
  Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme
  Buck, Matthias; Bouguet-Bonnet, Sabine; Pastor, Richard W.; MacKerell, Alexander D., Jr.
  Biophysical Journal (2006), 90 (4), L36-L38CODEN: BIOJAU; ISSN:0006-3495. (Biophysical Society)
  The recently developed CMAP correction to the CHARMM22 force field (C22) is evaluated from 25 ns mol. dynamics simulations on hen lysozyme. Substantial deviations from exptl. backbone root mean-square fluctuations and N-H NMR order parameters obtained in the C22 trajectories (esp. in the loops) are eliminated by the CMAP correction. Thus, the C22/CMAP force field yields improved dynamical and structural properties of proteins in mol. dynamics simulations.
59. 59
  Durell, S. R.; Brooks, B. R.; Ben-Naim, A. Solvent-induced forces between two hydrophilic groups. J. Phys. Chem. 1994, 98, 2198– 2202, DOI: 10.1021/j100059a038
  59
  Solvent-Induced Forces between Two Hydrophilic Groups
  Durell, Stewart R.; Brooks, Bernard R.; Ben-Naim, Arieh
  Journal of Physical Chemistry (1994), 98 (8), 2198-202CODEN: JPCHAX; ISSN:0022-3654.
  Mol. dynamics simulations were used to calc. the force between two simple hydrophilic solutes in dil. aq. soln. The "solutes" were two water mols. in the same relative orientation as the next-nearest neighbors in hexagonal ice I. Both the direct and solvent-induced contributions to the force were calcd. as a function of sepn. distance. The total force between the solutes was found to be most attractive at 5.0 Å (-1.6 kcal/mol/Å). The potential of mean force had a min. at 4.3 Å, which is 0.2 Å closer than the next-nearest-neighbor distance in ice. A parallel set of simulations were conducted with the partial charges on the "solutes" removed to examine hydrophobic analogs. In this case, the total force was most attractive at 3.5 Å (-0.9 kcal/mol/Å), and the min. of the potential was at the contact distance of 3.2 Å. In agreement with earlier predictions, the max. solvent-induced contribution to the potential was ca. 4 times more neg. for the hydrophilic "solutes" than for the hydrophobic ones. These differences are shown to be due predominantly to a solvent water mol. which simultaneously hydrogen bonds to both hydrophilic "solutes". The results support earlier assertions that solvent-induced interactions between polar amino acid residues are more important in protein folding and stability than generally considered.
60. 60
  Neria, E.; Fischer, S.; Karplus, M. Simulation of activation free energies in molecular systems. J. Chem. Phys. 1996, 105, 1902– 1921, DOI: 10.1063/1.472061
  60
  Simulation of activation free energies in molecular systems
  Neria, Eyal; Fischer, Stefan; Karplus, Martin
  Journal of Chemical Physics (1996), 105 (5), 1902-1921CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  A method is presented for detg. activation free energies in complex mol. systems. The method relies on knowledge of the min. energy path and bases the activation free energy calcn. on moving along this path from a min. to a saddle point. Use is made of a local reaction coordinate which describes the advance of the reaction in each segment of the min. energy path. The activation free energy is formulated as a sum of two terms. The first is due to the change in the local reaction coordinate between the endpoints of each segment of the path. The second is due to the change in direction of the min. energy path between consecutive segments. Both contributions can be obtained by mol. dynamics simulations with a constraint on the local reaction coordinate. The method is illustrated by applying it to a model potential and to the C7eq to C7ax transition in the alanine dipeptide. It is found that the term due to the change of direction in the reaction path can make a substantial contribution to the activation free energy.
61. 61
  Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable water model for the coarse-grained MARTINI force field. PLoS computational biology 2010, 6, e1000810, DOI: 10.1371/journal.pcbi.1000810
  61
  Polarizable water model for the coarse-grained MARTINI force field
  Yesylevskyy Semen O; Schafer Lars V; Sengupta Durba; Marrink Siewert J
  PLoS computational biology (2010), 6 (6), e1000810 ISSN:.
  Coarse-grained (CG) simulations have become an essential tool to study a large variety of biomolecular processes, exploring temporal and spatial scales inaccessible to traditional models of atomistic resolution. One of the major simplifications of CG models is the representation of the solvent, which is either implicit or modeled explicitly as a van der Waals particle. The effect of polarization, and thus a proper screening of interactions depending on the local environment, is absent. Given the important role of water as a ubiquitous solvent in biological systems, its treatment is crucial to the properties derived from simulation studies. Here, we parameterize a polarizable coarse-grained water model to be used in combination with the CG MARTINI force field. Using a three-bead model to represent four water molecules, we show that the orientational polarizability of real water can be effectively accounted for. This has the consequence that the dielectric screening of bulk water is reproduced. At the same time, we parameterized our new water model such that bulk water density and oil/water partitioning data remain at the same level of accuracy as for the standard MARTINI force field. We apply the new model to two cases for which current CG force fields are inadequate. First, we address the transport of ions across a lipid membrane. The computed potential of mean force shows that the ions now naturally feel the change in dielectric medium when moving from the high dielectric aqueous phase toward the low dielectric membrane interior. In the second application we consider the electroporation process of both an oil slab and a lipid bilayer. The electrostatic field drives the formation of water filled pores in both cases, following a similar mechanism as seen with atomistically detailed models.
62. 62
  Berendsen, H. J.; Postma, J. P.; van Gunsteren, W. F.; Hermans, J. Intermolecular forces; Springer, 1981; pp 331– 342.
  There is no corresponding record for this reference.
63. 63
  De Jong, D. H.; Baoukina, S.; Ingólfsson, H. I.; Marrink, S. J. Martini straight: Boosting performance using a shorter cutoff and GPUs. Comput. Phys. Commun. 2016, 199, 1– 7, DOI: 10.1016/j.cpc.2015.09.014
  63
  Martini straight: Boosting performance using a shorter cutoff and GPUs
  de Jong, Djurre H.; Baoukina, Svetlana; Ingolfsson, Helgi I.; Marrink, Siewert J.
  Computer Physics Communications (2016), 199 (), 1-7CODEN: CPHCBZ; ISSN:0010-4655. (Elsevier B.V.)
  In mol. dynamics simulations, sufficient sampling is of key importance and a continuous challenge in the field. The coarse grain Martini force field has been widely used to enhance sampling. In its original implementation, this force field applied a shifted Lennard-Jones potential for the non-bonded van der Waals interactions, to avoid problems related to a relatively short cutoff. Here we investigate the use of a straight cutoff Lennard-Jones potential with potential modifiers. Together with a Verlet neighbor search algorithm, the modified potential allows the use of GPUs to accelerate the computations in Gromacs. We find that this alternative potential has little influence on most of the properties studied, including partitioning free energies, bulk liq. properties and bilayer properties. At the same time, energy conservation is kept within reasonable bounds. We conclude that the newly proposed straight cutoff approach is a viable alternative to the std. shifted potentials used in Martini, offering significant speedup even in the absence of GPUs.
64. 64
  Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101, DOI: 10.1063/1.2408420
  64
  Canonical sampling through velocity rescaling
  Bussi, Giovanni; Donadio, Davide; Parrinello, Michele
  Journal of Chemical Physics (2007), 126 (1), 014101/1-014101/7CODEN: JCPSA6; ISSN:0021-9606. (American Institute of Physics)
  The authors present a new mol. dynamics algorithm for sampling the canonical distribution. In this approach the velocities of all the particles are rescaled by a properly chosen random factor. The algorithm is formally justified and it is shown that, in spite of its stochastic nature, a quantity can still be defined that remains const. during the evolution. In numerical applications this quantity can be used to measure the accuracy of the sampling. The authors illustrate the properties of this new method on Lennard-Jones and TIP4P water models in the solid and liq. phases. Its performance is excellent and largely independent of the thermostat parameter also with regard to the dynamic properties.
65. 65
  Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182– 7190, DOI: 10.1063/1.328693
  65
  Polymorphic transitions in single crystals: A new molecular dynamics method
  Parrinello, M.; Rahman, A.
  Journal of Applied Physics (1981), 52 (12), 7182-90CODEN: JAPIAU; ISSN:0021-8979.
  A Lagrangian formulation is introduced; it can be used to make mol. dynamics (MD) calcns. on systems under the most general, externally applied, conditions of stress. In this formulation the MD cell shape and size can change according to dynamic equations given by this Lagrangian. This MD technique was used to the study of structural transitions of a Ni single crystal under uniform uniaxial compressive and tensile loads. Some results regarding the stress-strain relation obtained by static calcns. are invalid at finite temp. Under compressive loading, the model of Ni shows a bifurcation in its stress-strain relation; this bifurcation provides a link in configuration space between cubic and hexagonal close packing. Such a transition could perhaps be obsd. exptl. under extreme conditions of shock.
66. 66
  Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: a linear constraint solver for molecular simulations. Journal of computational chemistry 1997, 18, 1463– 1472, DOI: 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H
  66
  LINCS: a linear constraint solver for molecular simulations
  Hess, Berk; Bekker, Henk; Berendsen, Herman J. C.; Fraaije, Johannes G. E. M.
  Journal of Computational Chemistry (1997), 18 (12), 1463-1472CODEN: JCCHDD; ISSN:0192-8651. (Wiley)
  We present a new LINear Constraint Solver (LINCS) for mol. simulations with bond constraints using the enzyme lysozyme and a 32-residue peptide as test systems. The algorithm is inherently stable, as the constraints themselves are reset instead of derivs. of the constraints, thereby eliminating drift. Although the derivation of the algorithm is presented in terms of matrixes, no matrix matrix multiplications are needed and only the nonzero matrix elements have to be stored, making the method useful for very large mols. At the same accuracy, the LINCS algorithm is 3-4 times faster than the SHAKE algorithm. Parallelization of the algorithm is straightforward.
67. 67
  Miyamoto, S.; Kollman, P. A. Settle: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. Journal of computational chemistry 1992, 13, 952– 962, DOI: 10.1002/jcc.540130805
  67
  SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models
  Miyamoto, Shuichi; Kollman, Peter A.
  Journal of Computational Chemistry (1992), 13 (8), 952-62CODEN: JCCHDD; ISSN:0192-8651.
  An anal. algorithm, called SETTLE, for resetting the positions and velocities to satisfy the holonomic constraints on the rigid water model is presented. This method is based on the Cartesian coordinate system and can be used in place of SHAKE and RATTLE. The authors implemented this algorithm in the SPASMS package of mol. mechanics and dynamics. Several series of mol. dynamics simulations were carried out to examine the performance of the new algorithm in comparison with the original RATTLE method. SETTLE is of higher accuracy and is faster than RATTLE with reasonable tolerances by three to nine times on a scalar machine. The performance improvement ranged from factors of 26 to 98 on a vector machine since the method presented is not iterative.
68. 68
  Hub, J. S.; de Groot, B. L.; Grubmüller, H.; Groenhof, G. Quantifying artifacts in Ewald simulations of inhomogeneous systems with a net charge. J. Chem. Theory Comput. 2014, 10, 381– 390, DOI: 10.1021/ct400626b
  68
  Quantifying Artifacts in Ewald Simulations of Inhomogeneous Systems with a Net Charge
  Hub, Jochen S.; de Groot, Bert L.; Grubmueller, Helmut; Groenhof, Gerrit
  Journal of Chemical Theory and Computation (2014), 10 (1), 381-390CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  Ewald summation, which has become the de facto std. for computing electrostatic interactions in biomol. simulations, formally requires that the simulation box is neutral. For non-neutral systems, the Ewald algorithm implicitly introduces a uniform background charge distribution that effectively neutralizes the simulation box. Because a uniform distribution of counter charges typically deviates from the spatial distribution of counterions in real systems, artifacts may arise, in particular in systems with an inhomogeneous dielec. const. Here, we derive an anal. expression for the effect of using an implicit background charge instead of explicit counterions, on the chem. potential of ions in heterogeneous systems, which (i) provides a quant. criterion for deciding if the background charge offers an acceptable trade-off between artifacts arising from sampling problems and artifacts arising from the homogeneous background charge distribution, and (ii) can be used to correct this artifact in certain cases. Our model quantifies the artifact in terms of the difference in charge d. between the non-neutral system with a uniform neutralizing background charge and the real neutral system with a phys. correct distribution of explicit counterions. We show that for inhomogeneous systems, such as proteins and membranes in water, the artifact manifests itself by an overstabilization of ions inside the lower dielec. by tens to even hundreds kilojoules per mol. We have tested the accuracy of our model in mol. dynamics simulations and found that the error in the calcd. free energy for moving a test charge from water into hexadecane, at different net charges of the system and different simulation box sizes, is correctly predicted by the model. The calcns. further confirm that the incorrect distribution of counter charges in the simulation box is solely responsible for the errors in the PMFs.
69. 69
  Thurlkill, R. L.; Grimsley, G. R.; Scholtz, J. M.; Pace, C. N. pK values of the ionizable groups of proteins. Protein science 2006, 15, 1214– 1218, DOI: 10.1110/ps.051840806
  69
  pK values of the ionizable groups of proteins
  Thurlkill, Richard L.; Grimsley, Gerald R.; Scholtz, J. Martin; Pace, C. Nick
  Protein Science (2006), 15 (5), 1214-1218CODEN: PRCIEI; ISSN:0961-8368. (Cold Spring Harbor Laboratory Press)
  We have used potentiometric titrns. to measure the pK values of the ionizable groups of proteins in alanine pentapeptides with appropriately blocked termini. These pentapeptides provide an improved model for the pK values of the ionizable groups in proteins. Our pK values detd. in 0.1 M KCl at 25°C are: 3.67 ± 0.03 (α-carboxyl), 3.67 ± 0.04 (Asp), 4.25 ± 0.05 (Glu), 6.54 ± 0.04 (His), 8.00 ± 0.03 (α-amino), 8.55 ± 0.03 (Cys), 9.84 ± 0.11 (Tyr), and 10.40 ± 0.08 (Lys). The pK values of some groups differ from the Nozaki and Tanford (N&T) pK values often used in the literature: Asp (3.67 this work vs. 4.0 N&T); His (6.54 this work vs. 6.3 N&T); α-amino (8.00 this work vs. 7.5 N&T); Cys (8.55 this work vs. 9.5 N&T); and Tyr (9.84 this work vs. 9.6 N&T). Our pK values will be useful to those who study pK perturbations in folded and unfolded proteins, and to those who use theory to gain a better understanding of the factors that det. the pK values of the ionizable groups of proteins.
70. 70
  Tanokura, M. 1H-NMR study on the tautomerism of the imidazole ring of histidine residues: I. Microscopic pK values and molar ratios of tautomers in histidine-containing peptides. Biochimica et Biophysica Acta (BBA)-Protein Structure and Molecular Enzymology 1983, 742, 576– 585, DOI: 10.1016/0167-4838(83)90276-5
  There is no corresponding record for this reference.
71. 71
  Chiang, C.-M.; Chien, K.-Y.; Lin, H.-j.; Lin, J.-F.; Yeh, H.-C.; Ho, P.-l.; Wu, W.-g. Conformational change and inactivation of membrane phospholipid-related activity of cardiotoxin V from Taiwan cobra venom at acidic pH. Biochemistry 1996, 35, 9167– 9176, DOI: 10.1021/bi952823k
  71
  Conformational Change and Inactivation of Membrane Phospholipid-Related Activity of Cardiotoxin V from Taiwan Cobra Venom at Acidic pH
  Chiang, Chien-Min; Chien, Kun-Yi; Lin, Hai-jui; Lin, Ji-Fu; Yeh, Hsien-Chi; Ho, Pei-li; Wu, Wen-guey
  Biochemistry (1996), 35 (28), 9167-9176CODEN: BICHAW; ISSN:0006-2960. (American Chemical Society)
  The phospholipid binding activity of cardiotoxin V from Naja naja atra (CTX A5) was studied by use of Langmuir monolayers and found to exhibit pH-dependence in binding to phosphatidylcholine membrane with an apparent pKa around 6.0. Proton NMR investigation of the CTX A5 mol. in the presence of phosphatidylcholine micelles reveals a decrease in assocn. of CTX A5 with membranes at low pH as a result of the protonation of His-4 near the membrane binding site of loop I region of CTX. The pH-dependent binding can be attributed mainly, but not solely, to the change in charge content of the CTX mol. upon His-4 protonation at the membrane/water interface. This is shown by analyzing the pH- and ionic strength dependence of binding of CTXs to phospholipid monolayers according to Gouy-Chapman theory. The protonation of the His-4 residue also results in a local conformational change in the loop I region since the chem. shifts of amide protons for the amino acid residues from Cys-3 to Thr-14 are all found to vary as a function of pH with an apparent pKa similar to that of His-4. Interestingly, the effect is relayed to other amino acid residues in the structural core of the protein such as those in C-terminal (Lys-60, Cys-61, and Asn-62) and triple-stranded antiparallel β-sheet (Cys-22, Lys-24, Ala-25, Arg-38, and Ala-41) regions. An addnl. local conformational change in the mol. results around pH 5 as evidenced by CD spectroscopic studies, although this change does not affect the characteristic β-sheet and three-finger loop structure of CTX mol. as revealed by two-dimensional NOESY 1H NMR study. The latter conformational change at acidic pH, however, completely inactivates CTX-induced aggregation/fusion activity of sphingomyelin vesicles. The results suggest that deciphering the functional sites of CTXs on the basis of structure and dynamics detd. at low pH should be done with caution. Since 19 out of 44 CTX homologs with known amino acid sequence contain His-4, the effect of His-4 on the structure and function of CTX mols. is important and is discussed in terms of the diverse membrane targets of CTX subtypes. Also discussed is the pH-induced activation of snake venom proteins in the victim.
72. 72
  Chiang, C.-M.; Chang, S.-L.; Lin, H.-j.; Wu, W.-g. The role of acidic amino acid residues in the structural stability of snake cardiotoxins. Biochemistry 1996, 35, 9177– 9186, DOI: 10.1021/bi960077t
  72
  The Role of Acidic Amino Acid Residues in the Structural Stability of Snake Cardiotoxins
  Chiang, Chien-Min; Chang, Shou-Lin; Lin, Hai-jui; Wu, Wen-guey
  Biochemistry (1996), 35 (28), 9177-9186CODEN: BICHAW; ISSN:0006-2960. (American Chemical Society)
  PH-induced structural transitions occurred in cardiotoxin V from Naja naja atra (CTX A5) at two pH values as judged by the CD ellipticity around 195 nm: an increase in the β-sheet content occurred around pH 4 and followed by a decrease, therein, around pH 2. The pKa of three acidic amino acid residues in CTX A5, i.e., Glu-17, Asp-42, and Asp-59, were detd. to be 4.0, 3.2, and below 2.3, resp., by NMR spectroscopy. The low pKa value of Asp-59 implies salt bridge formation between Lys-2 and Asp-59. Thus, electrostatic interaction may stabilize the three loop structure in addn. to the hydrogen bonds between N- and C-termini of CTX mol. Second, 2,2,2-trifluoroethanol (TFE) and guanidinium chloride (GdmHCl) were found to induce α-helical and random coil formation, resp., in CTX A5 and eight other β-sheet CTXs. Comparison of the relative potencies of TFE and GdmHCl to induce structural changes suggests that the amino acid residue located at position 17 plays a role in the structural stability. Specifically, CTXs contg. neg. charged Glu-17 are least stable. It is suggested that Glu-17 may perturb the interaction between Lys-2 and Asp-59, and thus the overall stability of β-sheet, in the presence of denaturing reagent. In conclusion, the perturbed structural stability of CTXs may partially explain the lower activity CTX exhibits at acidic pH. A structural model to account for the unfolding and refolding of CTX mols. without the breaking of disulfide bonds is also proposed.
73. 73
  Webb, H.; Tynan-Connolly, B. M.; Lee, G. M.; Farrell, D.; O’Meara, F.; Søndergaard, C. R.; Teilum, K.; Hewage, C.; McIntosh, L. P.; Nielsen, J. E. Remeasuring HEWL pKa values by NMR spectroscopy: Methods, analysis, accuracy, and implications for theoretical pKa calculations. Proteins: Struct., Funct., Bioinf. 2011, 79, 685– 702, DOI: 10.1002/prot.22886
  73
  Remeasuring HEWL pKa values by NMR spectroscopy: methods, analysis, accuracy, and implications for theoretical pKa calculations
  Webb, Helen; Tynan-Connolly, Barbara Mary; Lee, Gregory M.; Farrell, Damien; O'Meara, Fergal; Sondergaard, Chresten R.; Teilum, Kaare; Hewage, Chandralal; McIntosh, Lawrence P.; Nielsen, Jens Erik
  Proteins: Structure, Function, and Bioinformatics (2011), 79 (3), 685-702CODEN: PSFBAF ISSN:. (Wiley-Liss, Inc.)
  Site-specific pKa values measured by NMR spectroscopy provide essential information on protein electrostatics, the pH-dependence of protein structure, dynamics and function, and constitute an important benchmark for protein pKa calcn. algorithms. Titrn. curves can be measured by tracking the NMR chem. shifts of several reporter nuclei vs. sample pH. However, careful anal. of these curves is needed to ext. residue-specific pKa values since pH-dependent chem. shift changes can arise from many sources, including through-bond inductive effects, through-space elec. field effects, and conformational changes. We have re-measured titrn. curves for all carboxylates and His 15 in Hen Egg White Lysozyme (HEWL) by recording the pH-dependent chem. shifts of all backbone amide nitrogens and protons, Asp/Glu side chain protons and carboxyl carbons, and imidazole protonated carbons and protons in this protein. We extd. pKa values from the resulting titrn. curves using std. fitting methods, and compared these values to each other, and with those measured previously by 1H NMR. This anal. gives insights into the true accuracy assocd. with exptl. measured pKa values. We find that apparent pKa values frequently differ by 0.5-1.0 units depending upon the nuclei monitored, and that larger differences occasionally can be obsd. The variation in measured pKa values, which reflects the difficulty in fitting and assigning pH-dependent chem. shifts to specific ionization equil., has significant implications for the exptl. procedures used for measuring protein pKa values, for the benchmarking of protein pKa calcn. algorithms, and for the understanding of protein electro-statics in general.
74. 74
  Lee, J.; Miller, B. T.; Damjanovic, A.; Brooks, B. R. Enhancing constant-pH simulation in explicit solvent with a two-dimensional replica exchange method. J. Chem. Theory Comput. 2015, 11, 2560– 2574, DOI: 10.1021/ct501101f
  74
  Enhancing Constant-pH Simulation in Explicit Solvent with a Two-Dimensional Replica Exchange Method
  Lee, Juyong; Miller, Benjamin T.; Damjanovic, Ana; Brooks, Bernard R.
  Journal of Chemical Theory and Computation (2015), 11 (6), 2560-2574CODEN: JCTCCE; ISSN:1549-9618. (American Chemical Society)
  We present a new method for enhanced sampling for const.-pH simulations in explicit water based on a two-dimensional (2D) replica exchange scheme. The new method is a significant extension of our previously developed const.-pH simulation method, which is based on enveloping distribution sampling (EDS) coupled with a one-dimensional (1D) Hamiltonian exchange method (HREM). EDS constructs a hybrid Hamiltonian from multiple discrete end state Hamiltonians that, in this case, represent different protonation states of the system. The ruggedness and heights of the hybrid Hamiltonian's energy barriers can be tuned by the smoothness parameter. Within the context of the 1D EDS-HREM method, exchanges are performed between replicas with different smoothness parameters, allowing frequent protonation-state transitions and sampling of conformations that are favored by the end-state Hamiltonians. In this work, the 1D method is extended to 2D with an addnl. dimension, external pH. Within the context of the 2D method (2D EDS-HREM), exchanges are performed on a lattice of Hamiltonians with different pH conditions and smoothness parameters. We demonstrate that both the 1D and 2D methods exactly reproduce the thermodn. properties of the semigrand canonical (SGC) ensemble of a system at a given pH. We have tested our new 2D method on aspartic acid, glutamic acid, lysine, a four residue peptide (sequence KAAE), and snake cardiotoxin. In all cases, the 2D method converges faster and without loss of precision; the only limitation is a loss of flexibility in how CPU time is employed. The results for snake cardiotoxin demonstrate that the 2D method enhances protonation-state transitions, samples a wider conformational space with the same amt. of computational resources, and converges significantly faster overall than the original 1D method.
75. 75
  Pezeshkian, W.; König, M.; Wassenaar, T. A.; Marrink, S. J. Backmapping triangulated surfaces to coarse-grained membrane models. Nat. Commun. 2020, 11, 2296, DOI: 10.1038/s41467-020-16094-y
  75
  Backmapping triangulated surfaces to coarse-grained membrane models
  Pezeshkian, Weria; Koenig, Melanie; Wassenaar, Tsjerk A.; Marrink, Siewert J.
  Nature Communications (2020), 11 (1), 2296CODEN: NCAOBW; ISSN:2041-1723. (Nature Research)
  Many biol. processes involve large-scale changes in membrane shape. Computer simulations of these processes are challenging since they occur across a wide range of spatiotemporal scales that cannot be investigated in full by any single current simulation technique. A potential soln. is to combine different levels of resoln. through a multiscale scheme. Here, we present a multiscale algorithm that backmaps a continuum membrane model represented as a dynamically triangulated surface (DTS) to its corresponding mol. model based on the coarse-grained (CG) Martini force field. Thus, we can use DTS simulations to equilibrate slow large-scale membrane conformational changes and then explore the local properties at CG resoln. We demonstrate the power of our method by backmapping a vesicular bud induced by binding of Shiga toxin and by transforming the membranes of an entire mitochondrion to near-at. resoln. Our approach opens the way to whole cell simulations at mol. detail.
76. 76
  Wolfram Research, Inc. Mathematica, Version 11.3. Champaign: IL, 2018.
  There is no corresponding record for this reference.
Supporting Information
Supporting Information
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.2c00516.
- The input files and parameters used for MD simulations in the presented work (ZIP)
- A Mathematica notebook with instructions and routines to fit V^MM (ZIP)
- Description of λ-potentials, the effect of neglecting the interpolation of Lennard-Jones interactions, titration results for Asp and Glu within the single-site representation, a comparison of pK_a values for HEWL obtained with various λ-dynamics-based constant pH methods, demonstration that charge interpolation requires a single evaluation of the electrostatic potential for both single- and multisite representations (PDF)
Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system: http://pubs.acs.org/page/copyright/permissions.html.

Scalable Constant pH Molecular Dynamics in GROMACSClick to copy article linkArticle link copied!

Journal of Chemical Theory and Computation

License Summary*

Abstract

License Summary*

License Summary*

License Summary*

Introduction

Theory

λ-Dynamics-Based Constant pH MD Simulations

λ-Dependent Potential Energy Function

Figure 1

Linear Interpolation of Potential Energy Functions

Linear Interpolation of Partial Charges

Multisite Representation of Chemically Coupled Titratable Sites

Figure 2

Methods

Simulated Systems

Constant pH MD Simulation Setups

Reference States and Force Field Correction Potentials

Buffer Particles

Analysis of the Constant pH Trajectories

Results and Discussion

Titration of Single Amino Acids

Figure 3

Titration of Proteins

Figure 4

Figure 5

Efficiency of the Implementation

Figure 6

Conclusions

Supporting Information

Terms & Conditions

Author Information

Acknowledgments

Appendix: Constraint Algorithm

References

Cited By

Journal of Chemical Theory and Computation

License Summary*

Article Views

Altmetric

Citations

Recommended Articles

Abstract

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

References

Supporting Information

Supporting Information

Terms & Conditions

Scalable Constant pH Molecular Dynamics in GROMACS
Click to copy article linkArticle link copied!