A pH Replica Exchange scheme in the Stochastic Titration constant

Subscriber access provided by OCCIDENTAL COLL

Molecular Mechanics

A pH Replica Exchange scheme in the Stochastic Titration constant-pH MD method Diogo Vila-Viçosa, Pedro B. P. S. Reis, António M. Baptista, Chris Oostenbrink, and Miguel Machuqueiro J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 25 Mar 2019 Downloaded from http://pubs.acs.org on March 27, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

A pH Replica Exchange scheme in the Stochastic Titration constant-pH MD method Diogo Vila-Vi¸cosa,∗,†,‡ Pedro B. P. S. Reis,†,‡ António M. Baptista,¶ Chris Oostenbrink,§ and Miguel Machuqueiro†,‡ †Centro de Qu´ımica e Bioqu´ımica, Departamento de Qu´ımica e Bioqu´ımica, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal ‡University of Lisboa, Faculty of Sciences, BioISI - Biosystems & Integrative Sciences Institute, Lisboa, Portugal ¶Instituto de Tecnologia Qu´ımica e Biológica António Xavier Universidade Nova de Lisboa, Av. da Rep´ ublica, 2780-157 Oeiras, Portugal §Department of Material Sciences and Process Engineering, Institute of Molecular Modeling and Simulation, University of Natural Resources and Life Sciences Vienna, Muthgasse 18, A-1190 Vienna, Austria E-mail: [email protected] Phone: +351-21-7500112. Fax: +351-21-7500088 Abstract Solution pH is a physicochemical property that has a key role in cellular regulation and its impact at the molecular level is often difficult to study by experimental methods. In this context, several theoretical methods were developed to study pH effects in macromolecules. The stochastic titration constant-pH molecular dynamics method (CpHMD) was developed by coupling molecular sampling methods, that are appropriate to study the conformational ensemble of biomolecules, with continuum electrostatics

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

approaches, that properly describe pH-dependent protonation states. However, in difficult cases, the protonation sampling can be too slow for the commonly accessible computational times. In this work, we combined a pH replica exchange scheme with this CpHMD method, exploring several optimization strategies and possible limitations.

1

Introduction

Over the last decades, several theoretical methods were developed to study the pH effects in macromolecules at the microscopic level. Continuum electrostatic methods are used to estimate the electrostatic potential generated by a solute embedded in a solution that is described implicitly by a dielectric constant. With this estimation, it is possible to calculate the electrostatic energy associated with a (de)protonation event and, consequently, a free energy of protonation (pKa ). This general approach can be implemented using simplified methods such as Gouy–Chapman (an analytical solution of the Poisson–Boltzmann equation) 1 or generalized Born model. 2 A more detailed estimation of protonation free energies can be obtained using the Poisson–Boltzmann (PB) equation. 3–10 These energy terms can then be used to sample protonation states using a Monte Carlo simulation. 11–13 The combination of these two methods allows for an accurate estimation of pKa values of titrable amino acid residues based on experimental (rigid) structures. 14 To obtain information on the protonation/conformation behavior it is also necessary to sample conformations using a simulation method such as molecular dynamics (MD) or Monte Carlo (MC). By setting a typical protonation state of a given pH value, one can sample protonation-dependent conformational spaces. However, this state is often unknown and not unique. Another possibility is to run MD simulations at all possible protonation states which is also a rather limited approach since the number of protonation states varies with 2N (N being the number of titrable groups) meaning that the number of combinations becomes too large for a relatively small number of titrable groups. A possible solution, is to take into account (de)protonation events and only sample the more relevant states. This idea inspired the so-called constant-pH 2


Page 2 of 31

Page 3 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MD methods (CpHMD) 15–30 in which the protonation states of titrable sites are dynamic, fluctuating according to the external pH and to the solute conformation as a discrete or continuous variable. The stochastic titration CpHMD method, developed by Baptista and coworkers 18,23 is based on a cyclic stop-and-go algorithm where the conformations are sampled with MD using a classical force field and protonation states are obtained from MC using PB-derived free-energy terms. Simulations of biological systems with standard MD simulations are often affected by sampling limitations (e. g. events that occur at too large timescales). Enhanced sampling methods coupled with MD simulations are often used to overcome these problems, by either still sampling from the proper ensemble, or allowing for its reconstruction a posteriori. In this context, several methods have appeared whose review is beyond the scope of this text (please see 31,32 and references therein for a detailed overview and discussion of enhanced sampling methods). CpHMD methods face similar hurdles related to the sampling efficiency of protonation states 26–29,33 and, therefore, also benefit from enhanced sampling techniques. In fact, the combination of accelerated MD with a CpHMD method was already proposed. 34 Also, the sampling of such simulations can be increased using a coarse grained model that allows to reach longer simulation times. 35 However, as mentioned, most enhanced sampling methods were developed for conformations and need some adjustments to deal with protonation states. In CpHMD simulations, this problem has been addressed using replica exchange techniques, originally developed for temperature where several replicas can exchange between different temperatures. 36 These exchanges allow for the mixing of more diverse conformational states that are more easily sampled at higher temperatures. This concept can also be generalized to exchange replicas between different Hamiltonians in single or multidimensional replica exchange schemes. 37 In the context of CpHMD simulations, replicas are allowed to switch between pH values mixing conformations (and protonation states) from different pH values. 26–29,33 The criterion to swap between pH values must be carefully determined in order to sample the proper pH dependent ensembles (see Theory

3



Page 4 of 31

and Methods). Using the enveloping distribution sampling (EDS) scheme, a Hamiltonian replica exchange method (HREM) was proposed. 38 In this case, replicas exchange between Hamiltonians corresponding to different EDS potential energy surfaces. 38 Furthermore, a 2D replica exchange method was developed, which combines exchanges between pH values and between Hamiltonians. 30 Here, we implement the pH replica exchange scheme (pHRE) proposed by Itoh et al. 26 in the context of the stochastic titration CpHMD method. 18,23 We test it using a model system (ethylenediamine - EDA) that has two experimentally known macroscopic pKa values 39,40 and a typical benchmark system for pKa calculations (hen-egg white lysozyme - HEWL). 41,42 We discuss the advantages and limitations of the pHRE enhanced sampling technique and test some of its parameters. As previously observed for other pHRE implementations 26,28 we were able to obtain identical titration curves with both methods (CpHMD and pHRE) and a faster protonation sampling for pHRE simulations of HEWL. We carefully investigate how the convergenve and sampling depends on the time between exchanges. As both the replica exchange method as well as the stochastic titration method have characteristic relaxation times between exchanges or protonation events, a particular focus is on the interplay between these times.

2 2.1

Theory and methods pH Replica Exchange

Let i and j be two replicas at pH values of pHl and pHm , respectively. In each replica exchange step the probability of accepting a trial exchange, Pacc , is evaluated between two states. Being a Metropolis Monte Carlo step, 43 this probability is given by:

Pacc

Pm (Ri )Pl (Rj ) = min 1, Pl (Ri )Pm (Rj )

4


(1)

Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where Pl is the equilibrium probability of state Ri characterized by a coordinate vector (qi ), a momentum vector (pi ), a protonation state (xi ), and a pH value (pHl ). For a system at constant temperature (T ), pressure (Pe) and pH, the equilibrium probability is given by:

Pl (Ri ) =

1 exp[−β(K(pi ) + U (qi , xi ) + kB T pHl N (xi ) ln 10 + PeV )] Ξl

(2)

where Ξl is the partition function, V is the volume, U is the potential energy, K is the kinetic energy (note that, in the stochastic titration constant-pH MD method, the kinetic energy does not depend on the protonation state since the titrable protons are always present either as dummy or fully charged atoms) and N (xi ) is the number of titrable groups that are protonated in state xi . For example, the exchange can correspond to moving the positions and momenta of replica i to pHm and of replica j to pHl . Combining the two equations above, Pacc is then given by:

Pacc = min{1, exp[−(pHm − pHl )(N (xi ) − N (xj )) ln 10]}

(3)

as previously reported. 26 Interestingly, we note that there are no energy dependencies and the acceptance probability only depends on the pH values of the two replicas and on the number of protonated groups. However, it is also possible to exchange both pH and protonation states from one replica to the other. That would mean that replica i (j) that is being simulated at pHl (pHm ) and has a N (xi ) (N (xj )) number of protons will be moved to pHm (pHl ) and will have N (xj ) (N (xi )) titrable protons. The acceptance probability is now given by:

Pacc = min{1, exp[β(−∆U (qi , xi → xj ) − ∆U (qj , xj → xi ))]}

(4)

where the explicit dependence on pH values and protonation states gets canceled. In this case, the acceptance probability depends on the potential energy difference of picking up a

5



Page 6 of 31

protonation state from a different conformation. In other words, it is an estimation of how well the two conformations can accommodate (in terms of potential energy) the protonation state of the each other. This potential energy term is, in general, very high which means that the probability of exchange could be rather low, rendering this exchange scheme very inefficient. There is one last possibility that corresponds to exchanging only the protonation state and keep the replicas at the same pH value. This acceptance probability is given by:

Pacc = min{1, exp[β(−∆U (qi , xi → xj )−∆U (qj , xj → xi ))−(pHm −pHl )(N (xi )−N (xj )) ln 10]} (5) where we have all previous dependencies combined: pH, number of protons and potential energy differences. However, as occurred in the second scenario, we still have the potential energy difference in our equation that hinders the efficiency of the replica exchange method. In this work, all pH replica exchange simulations were performed using equation 3 to calculate the acceptance probability. Being identical to the schemes used in References 26 and 28, we also note that, in this scenario, the acceptance probabilities will be higher and we could obtain a better replica mixing.

2.2

Constant-pH MD and pH Replica Exchange Settings

The stochastic titration constant-pH molecular dynamics (MD) method (CpHMD) allows for the treatment of the pH value as an external parameter in a MD simulation. 18,23 In this method, a molecular mechanics/molecular dynamics (MM/MD) simulation is interrupted at regular time intervals (τprt ) and new protonation states are obtained from a Monte Carlo (MC) simulation using Poisson–Boltzmann (PB) derived free-energy terms. 18,23 The implementation of the pH replica exchange method was achieved using a Python-based wrapper to the previous CpHMD implementation. 18,23 Both pHRE and CpHMD codes are available from the authors upon request. In this approach, it is necessary to define the

6


Page 7 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


frequency of the exchange attempts (in our case, we perform a replica exchange step every τRE ). In most of our simulations τprt is 20 ps and two values of τRE are tested for both systems: 20 and 100 ps. The only exception is in Section 3.3, where these two parameters are tested in the EDA system with values ranging from 2 to 100 ps. The exchanges are only attempted between neighboring pH values. An alternating exchange scheme was used: at odd exchange trials, exchanges between the (2i+1)th and (2(i+1))th pH values are attempted and at even exchange trials, exchanges between the (2i)th and the (2i+1)th pH values are attempted, with i being an auxiliary index to generate all possible pairs (see Figure S1 in Supporting Information for an example). This means that the effective time between two trials of the same pH values is 2τRE , avoiding immediate backswitches. 44 The ethylenediamine (EDA) simulations were performed at 12 equally-spaced pH values ranging from 6.0 to 11.5 for 200 ns (1 ns was discarded as equilibration) and both amino groups were allowed to titrate. Hen-egg white lysozyme (HEWL) simulations were performed for 100 ns, however, only the 5–30 ns segment was used for analysis to ensure structural integrity (see results and discussion). All acidic (Glu and Asp) and His residues and both termini were allowed to titrate. All simulations (CpHMD and pHRE) were performed in triplicate with 12 equally-spaced pH values ranging from 1.0 to 6.5. In both MM/MD and PB/MC steps the temperature was set to 310 K and the ionic strength to 0.1 M. All analyses were performed using GROMACS 4.0.7 45 and in-house software. Uncertainty values for single replicates were obtained using the autocorrelation function to compute the number of independent blocks. 46 For HEWL (where we use replicates), we report the standard error calculated from both the averages of each replicate and the autocorrelation function. For the pKa calculations of individual groups in HEWL, a leave one out strategy was used to calculate statistical uncertainty values similar to the reported in Ref. 47. All reported uncertainty values correspond to 2 × standard deviations.

7



2.3

Molecular mechanics/molecular dynamics

All molecular dynamics (MD) simulations were performed with an integration step of 2 fs using a modified version of GROMACS 4.0.7 45 that allows the user to set the ionic strength value. 23 The GROMOS 54a7 48,49 force field was used for HEWL, whereas the parameters for EDA were taken from reference 50. Water molecules were treated using the simple point charge (SPC) model. 51 Nonbonded interactions were treated with a twin-range cutoff scheme where all pairs with distance below 0.8 nm were updated every step and between 0.8 and 1.4 nm every 5 steps. Beyond 1.4 nm all van der Waals interactions were truncated and the electrostatics were complemented with a generalized reaction field 52 with a dielectric constant of 61. 53 In the NPT ensemble, temperature and pressure were kept constant using the v-rescale 54 thermostat (coupling constant of 0.1) and Parrinello-Rahman 55 barostat (coupling constant of 0.5 and isothermal compressibility of 4.5 × 10−5 bar−1 ). Both systems were solvated in a rhombic dodecahedral box with periodic boundary conditions with 579 water molecules for EDA and 5685 for HEWL.

2.4

Poisson–Boltzmann/Monte Carlo

Poisson–Boltzmann (PB) calculations were performed using DelPhi 5.1. 56 A cubic grid with 51 or 81 points (for EDA and HEWL, respectively) was used to ensure an initial grid space of ∼0.1 nm. This first step was always followed by a focusing procedure 57 on the group of interest (titrable group) using a grid with the same number of points but a grid spacing of ∼0.025 nm. In the convergence procedure, the threshold was set to 0.01 kB T /e. For the definition of the molecular surface, a probe with 0.14 nm and an ion exclusion layer of 0.2 nm were used. All charges were taken from the MM force field, whose Lennard-Jones parameters were also used to derive atomic radii. 14 The inner (outer) dielectric constant was set to 2 (80) in agreement with previous works. 18,23 The Monte Carlo (MC) sampling of protonation states was performed using PB-derived protonation free energies. These simulations were performed using PETIT 1.6.1 software 13 8


Page 8 of 31

Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with 105 cycles. Each cycle is defined as an attempt to change the protonation state of all titrable groups and pairs (only considered for interaction energies larger than 2 pK units).

3

Results and discussion

The pH replica exchange scheme (pHRE) here implemented within the stochastic titration constant-pH MD method (CpHMD) was tested in two systems: ethylenediamine (EDA) and hen-egg-white lysozyme (HEWL). The method efficiency was evaluated in terms of its ability to mix conformations from different pH values and the results were also compared with CpHMD simulations and experiments.

3.1

pHRE efficiency

The efficiency of a replica exchange scheme can be evaluated in several ways. Usually, it is important to have a high exchange rate (number of accepted exchanges per trial) and a significant number of roundtrips (a roundtrip is defined as two complete trips across all pH values - from the lower to the higher and vice versa). We determine the actual probability of exchange (Pexc ) that we observe in our simulations (Figure 1). In this measure, a Pexc of 1 means that all exchanges are accepted. Due to the alternating scheme of exchange trials, this represents 25 exchanges per ns for τRE = 20 ps and 5 for τRE = 100 ps (see Methods). As expected, Pexc does not vary with different τRE values, since Pacc is independent from the frequency of trial (equation 3). In our systems, we observe a higher Pexc value for EDA simulations in the entire simulated pH range. This is expected since Pacc depends on the difference of titrable protons between the conformations from two pH values (equation 3). Since the ∆pH was equal in the two systems, the one with less titrable sites (EDA) is expected to show higher Pexc values (∆N values will, on average, be smaller). This result shows that the pHRE scheme mixes conformations more efficiently in EDA simulations when compared with HEWL. However, if necessary, Pexc can be increased by decreasing ∆pH (equation 3).

9



1

20 ps 100 ps

Pexc

0.8 0.6 0.4 0.2

6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5

pH 1

20 ps 100 ps

0.8 Pexc

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

0.6 0.4 0.2

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5

pH

Figure 1: Probability of exchange (Pexc ) between all pairs of pH values in all pHRE simulations of both EDA (top) and HEWL (bottom). Each bar corresponds to the Pexc between two adjacent pH values.

10


Page 10 of 31

Page 11 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The differences in the obtained Pexc between pH values are also explained by its dependence on the difference in the number of protons. In fact, Pexc is strongly correlated with the slope of the titration curve (described below and depicted in Figure 2). For example, the EDA titration curve has almost a plateau between pH 8.0 and 9.0 which corresponds to a Pexc of almost 1 (when the difference in the number of titrable protons is 0, Pacc is 1). For HEWL, the lower values of Pexc were observed for the pH values where the slope of the titration curve is higher (∼ 3 − 4). When the shape of the titration curve is known a priori, which is often not the case, this information can be used to perform a better choice of the pH values. In fact, this choice could be optimized based on short or previous simulations to obtain a constant Pexc for all simulated pH range (see Figure S2 in Supporting Information for an example of a pHRE simulation with an approximately constant Pexc ). The number of roundtrips was also evaluated and, for EDA simulations, we observed 0.46 ± 0.01 roundtrips per ns for τRE = 20 ps and 0.13 ± 0.01 for τRE = 100 ps. This difference indicates that a lower τRE value leads to a more efficient mixing between pH values and, probably, to a better sampling. For HEWL simulations very few roundtrips were observed. Some of the replicas were simulated at all pH values but most of them were restricted to a range of ∼ 2 − 4 pH units. This represents a significant difference between the two systems: in EDA simulations all replicas are virtually indistinguishable in terms of sampled pH values whereas in HEWL there is a significant mixing between pH values, but not uniformly over all replicas (Figures S3–S10 in Supporting Information). The efficiency of this pHRE method is dependent on the system (on its number of titrable groups and their pKa values) and the choice of pH values to run the simulations. Due to the larger number of titrable protons in HEWL, the efficiency of the pHRE method is lower. However, since the mixing between replicas is significant and we are interested in all pH values, this lower efficiency is not as problematic as it can be for other replica exchange schemes. Also, a large Pexc and number of roundtrips can introduce artifacts if our system has no time to structurally adapt to the new pH value. In other words, to avoid artifacts the

11



effective time between exchanges should be larger than the relaxation timescales.

3.2

pHRE reproduces CpHMD titration curves

The total titration curves obtained from the pHRE method are able to properly reproduce the results from the stochastic titration CpHMD in both EDA and HEWL systems (Figure 2). For EDA, the two protonation transitions are well captured with both methods and all 2

Robinson 1959 Haynes CpHMD τRE 100 ps τRE 20 ps

Total charge

1.5

1

0.5

0 6

7

8

9

10

11

pH Tanford 1972 CpHMD τRE 100 ps τRE 20 ps

18 16 Total charge

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

14 12 10 8 1

2

3

4

5

6

pH

Figure 2: Titration curves obtained for EDA (top) and HEWL (bottom) with both CpHMD and pHRE methods.

curves are mostly identical (CpHMD and pHRE), with small differences at pH 6.0 and 6.5 (see Section 3.3). This is reflected in the similarity between the calculated pKa values 12


Page 12 of 31

Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(Table 1). The small differences observed when comparing with experimental values can Table 1: Experimental 39,40 and simulated pKa values of EDA. Experimental

pHRE

Ref. 39 Ref. 40 CpHMD

τRE = 100 ps

τRE = 20 ps

pKa 1

7.00

6.86

6.96

6.94

6.93

pKa 2

10.09

9.92

10.21

10.21

10.22

∆pKa

3.09

3.06

3.25

3.27

3.29

be attributed to force field limitations since the used parameters were not calibrated to describe (de)protonation events in EDA. Furthermore, for two strongly coupled titrable groups, there are significant electronic effects contributing to their pKa that are not easily described with point charges. Even though it is often possible to adjust charge sets to obtain the experimental values, these parameters would hardly be transferable to other systems. The HEWL calculated titration curves are also indistinguishable (the observed differences are within the error bars for all pH values). Moreover, as previously observed, 58,59 there is a generally good agreement with the experimental curve and a small deviation in the more acidic region (see below). The agreement between all computational estimations is also reflected in the pKa values of the individual residues (Table 2; titration curves in Figures S11– S14 in Supporting Information). The pKa of most residues is properly captured by both methodologies compared with experimental values 41 as previously observed for CpHMD simulations using GROMOS 43a1 and 53a6 force fields. 58,59 In fact, the estimations for the most challenging residues with CpHMD (Glu35, Asp48, and Asp66) are not improved by pHRE. Overall, the root-mean-square error (RMSE) values are ∼0.83 pKa units in agreement with Ref. 59 (0.79 pKa units for GROMOS 53A6 force field). The Glu35 residue pKa value is often difficult to predict using computational methodologies 14,58,59 due to its location, inside the hydrophobic active site of HEWL. A proper description of the equilibrium between buried and solvent exposed conformations of Glu35 is crucial to obtain a correct pKa shift. In CpHMD simulations this conformational equilibrium 13



Page 14 of 31

Table 2: Experimental 41,42 and simulated pKa values of all titrable groups in HEWL. The root mean square error (RMSE) values are calculated comparing our pKa values with both experimental works. For the C-Ter pKa value from Ref. 42, the lowest value is used since it is more similar to the one from Ref. 41, allowing for easier comparison with previous results using CpHMD. Residue Bartik et al. 41

Webb et al. 42

CpHMD

τRE = 100 ps

τRE = 20 ps

Glu7

2.85

2.6

3.38 ± 0.07

3.36 ± 0.08

3.34 ± 0.04

His15

5.36

5.5

5.49 ± 0.06

5.43 ± 0.08

5.43 ± 0.04

Asp18

2.66

2.8

3.64 ± 0.10

3.65 ± 0.05

3.57 ± 0.07

Glu35

6.20

6.1

5.46 ± 0.37

5.68 ± 0.13

5.52 ± 0.10

Asp48

1.60

1.4

2.09 ± 0.41

2.11 ± 0.12

1.95 ± 0.15

Asp52

3.68

3.6

3.92 ± 0.22

4.01 ± 0.08

3.85 ± 0.14

Asp66

0.90

1.2

3.12 ± 0.28

3.20 ± 0.05

3.15 ± 0.04

Asp87

2.07

2.2

2.30 ± 0.15

2.26 ± 0.05

2.31 ± 0.08

Asp101

4.09

4.5

3.90 ± 0.24

3.98 ± 0.08

3.77 ± 0.10

Asp119

3.20

3.5

2.72 ± 0.11

2.73 ± 0.10

2.80 ± 0.08

C-Ter RMSE

2.75 –

2.7 –

3.20 ± 0.09 0.83|0.81

3.22 ± 0.08 0.83|0.81

3.24 ± 0.08 0.82|0.80

14


Page 15 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


is poorly described and needs to be captured by running three independent simulations at each pH value (Figure 3). With pHRE, due to the replica mixing, the use of multiple runs becomes partially unnecessary, since the obtained pKa values show smaller dispersion. This smaller dispersion is more evident for a τRE value of 20 ps due to a higher exchange frequency. Consequently, the pKa uncertainty values in our calculations tend to increase with τRE (note that a CpHMD simulation is equivalent to τRE → ∞). All pKa calculations of titrable residues in HEWL were performed using short segments (30 ns) ensuring the protein structural integrity. Nevertheless, we extended all these simulations to 100 ns and evaluated the conformational deviations from the X-ray structure (Figure 4). After 30 ns, we observe a significant increase in RMSD for some replicas, triggered by extremely low pH values (1–2). In pHRE simulations, the altered conformations spread over multiple pH values due to the replica mixing protocol. The conformational transition captured by the RMSD increase corresponds to an opening of the catalytic cleft which might be correlated with the initial steps of the experimentally observed HEWL acidic unfolding. In fact, guanidinium denaturation studies show that the protonation of some acidic groups with lower pKa values is key for the destabilization of HEWL at very low pH. 60,61 In our simulations, we bring the system to visit very low pH values through the pHRE method (trial scheme and equation 3), allowing the protonation of all residues, which may be inducing the aforementioned low pH denaturation transition.

3.3

A small discrepancy in the EDA titration curve

As previously mentioned, most of our results are consistent between CpHMD and pHRE. However, in the EDA titration curve, there is a small discrepancy between the two methods. This small, yet significant, difference only occurs at the two lowest simulated pH values (6.0 and 6.5; Figure 2 and 5a). At low pH, where the EDA molecule is almost constantly doubly protonated, the trans-gauche equilibrium is difficult to attain since the gauche conformation is highly unstable. In pHRE simulations, due to a high acceptance probability (see above), 15



1.0

Occupation

0.8

CpHMD

0.6 0.4 pKa = 5.46 ± 0.37

0.2

1.0

Occupation

0.8

τRE 100 ps

0.6 0.4 pKa = 5.68 ± 0.13

0.2

1.0 0.8 Occupation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

τRE 20 ps

0.6 0.4 pKa = 5.52 ± 0.10

0.2 0.0 3.0

4.0

5.0 pH

6.0

Figure 3: Titration curves obtained for Glu35 residue in HEWL with both CpHMD and pHRE methods.

16


Page 16 of 31

1.6

RMSD / nm

1.2

0.8

0.4 Replicate 1 1.6

1.2 RMSD / nm

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.8

0.4 Replicate 2 1.6

1.2 RMSD / nm

Page 17 of 31

pH 1.0 pH 2.0 pH 3.0 pH 4.0 pH 5.0 pH 6.0

0.8

0.4

0

Replicate 3 40

60 time / ns

80

100

Figure 4: RMSD values over time of pHRE simulations (all 3 replicates with τRE = 20 ps). Between 0 and 30 ns all average RMSD values are below 0.6 (data not shown).

17



Ethylenediamine − τprt 20 ps (zoom)

Ethylenediamine − τprt 20 ps (τRE 20 ps)

2

2

b) Protonation

Total charge

a) 1.9

1.9

1.8

CpHMD pH 5.5−11.5 pH 5−6 pH 6−11.5

1.8 CpHMD τRE 20 ps 6

6.25 pH

6.5

5.5

Ethylenediamine − τprt 20 ps (τRE variation)

6 pH

2

d) Protonation

c) 1.9

1.8

6.5

Ethylenediamine − τprt 2 ps (τRE variation)

2

Protonation

1.9

CpHMD τRE 100 ps τRE 40 ps τRE 20 ps 6

1.8

6.25 pH

CpHMD τRE 20 ps τRE 10 ps τRE 2 ps

6.5

6

6.25 pH

6.5

Ethylenediamine − τprt 100 ps 2

e) Protonation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 31

1.9

1.8 CpHMD τRE 100 ps 6

6.25 pH

6.5

Figure 5: Zoom of the titration curves obtained for EDA with both CpHMD and pHRE methods. The effect of multiple combinations of pH values, τprt and τRE is tested: (a) our “standard” approach (τprt = τRE = 20 ps); (b) the choice of a different set of pH values; the τRE variation for τprt = 20 ps (c) and τprt = 2 ps (d); and (e) τprt = τRE = 100 ps.

18


Page 19 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


at lower pH values the gauche configurations are being oversampled by replicas coming from higher pH values. This effect still holds by adding an extra pH (5.5) but vanishes if only low values are used (5.0, 5.5, and 6.0; Figure 5b). The oversampling of gauche conformations at pH 6.0 is also mitigated if τRE is increased (pHRE step is performed less frequently; Figures 5c and d), which is expected since the limit case (for an infinite τRE value) corresponds to a CpHMD simulation. The observed discrepancies are not specific for τRE = 20 ps, since they hold for other τRE values, being more evident when τprt = τRE (Figures 5c-e). Interestingly, by increasing the τRE values, the average protonation values of CpHMD and pHRE seem to converge to a common value. In fact, due to a particularly strong protonation/conformation coupling in the EDA molecule, there is a significant dependence of the average protonation on the τprt value. Nevertheless, these differences are similar to our observations in the pHRE and CpHMD comparison, further supporting the pHRE reliability. Care should be taken to allow the system sufficient time to relax the conformation and protonation states between replica exchange switches.

3.4

pHRE method accelerates protonation convergence

To evaluate the convergence speedup in pHRE simulations, we calculated the average root mean square (rms) value of the number of protons between a single replicate and the ensemble averages over time (Figure 6). For EDA, all deviations converge to 0, since only one replicate was performed. In this case, the protonation values converge extremely fast with no significant differences between the two methods, which was expected for such a simple system. For HEWL, the observed deviations are larger for CpHMD simulations due to larger differences in average protonation values between replicates, as previously mentioned (Figure 3). Also, the error bars for pHRE, are significantly smaller even for short simulation timescales (∼ 8 ns). These observations support that the use of replicates is almost redundant in pHRE simulations since a single run is able to properly capture the protonation of HEWL titrable residues. In essence, according to our results for HEWL and in agreement with previous reports, 26–29 19



0.01

CpHMD τRE 100 ps τRE 20 ps

rms / # protons

0.008 0.006 0.004 0.002 0 0

40

80 120 time / ns

0.2

160

200

25

30

CpHMD τRE 100 ps τRE 20 ps

0.16 rms / # protons

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

0.12

0.08

5

10

15 20 time / ns

Figure 6: Average root mean square (rms) value of the number of protons between a single replicate and the ensemble average for EDA (top) and HEWL (bottom). Error bars for HEWL correspond to the standard error over replicates.

20


Page 21 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the real advantage of using pHRE coupled to the stochastic titration CpHMD method relies on the need for fewer replicates to sufficiently sample the protein conformation/protonation space.

4

Conclusions

Many enhanced sampling techniques are nowadays well established to speed up MD simulations of biological systems. However, in the context of protonation sampling, namely in constantpH MD simulations, these are still emerging techniques that need validation. Here, we present the implementation of a pH replica exchange technique within the framework of the stochastic titration CpHMD method. 18,23 Several approaches are discussed, namely the exchange of pH value, protonation state or both. The first option is, in principle, more efficient and was adopted here, as proposed in Reference 26. Our choice seems adequate to be used with this CpHMD method since we obtained identical titration curves for both a simple and a more complex system. Most differences between CpHMD and pHRE simulations are within the statistical uncertainty. The efficiency of the pHRE method depends fundamentally on two parameters: τRE (the inverse of the exchange attempts frequency) and simulated pH values. Our overall results suggest that τRE should be as low as possible, maximizing the number of transitions. However, this value can originate small artifacts in small systems with strong protonation/conformation coupling such as ethylenediamine. This effect vanishes with larger τRE value and is also absent in HEWL simulations. This means that τRE should be long enough for the system to relax towards the new pH-value. In systems with a massive number of conformation states, but no major slow conformational changes, we suggest to use τRE = 20 ps. Regarding the choice of the pH values, it should depend on the total number of the titration sites and on the slope of the titration curve. The optimization of the pH spacing, to obtain constant exchange probabilities, can be achieved through a fast mimicking approach 62 or by running short simulations and manually adjust the pH values.

21



Page 22 of 31

As previously observed, 26–29 the major advantage of this enhanced sampling technique is to sample the conformational/protonation space faster than standard methods. Here, we showed that, for long enough simulations of the two test systems, the sampled protonation space is identical for both methods. The main difference between the two approaches is that with pHRE we reach the final results in less simulation and computational time. The pH replica exchange scheme in the stochastic titration constant-pH MD method seems to provide accurate results faster than the standard CpHMD. Its reliability allows it to be tested in larger and more complex systems, including membrane environments. 63

Acknowledgement We thank Jan W. Perthold, Pascal T. Merz and Tomás F. D. Silva for fruitful discussions. This work was financially supported by: Project LISBOA-01-0145-FEDER-007660 (Microbiologia Molecular, Estrutural e Celular) funded by FEDER funds through COMPETE2020 - Programa Operacional Competitividade e Internacionalizao (POCI) and by national funds through FCT - Fundao para a Cincia e a Tecnologia; The authors also acknowledge financial support from Funda¸cão para a Ciência e a Tecnologia through grant SFRH/BPD/110491/2015, postdoctoral grant under project IF/00069/2014/CP1216/CT0006, and project grants UID/MULTI/00612/2013, UID/MULTI/00612/2019, UID/MULTI/04046/2019, PTDC/QUI-OUT/29441/2017, and PTDC/QEQCOM/5904/2014.

Dedication Dedicated to Prof. Maria José Calhorda on the occasion of her 70th birthday and retirement.

22


Page 23 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Supporting Information Available A schematic representation of the used exchange scheme. A graphical representation of the Pexc for the optimized pHRE simulation of EDA. Temporal variation of pH values sampled by each replica in all pHRE simulations. Titration curves for individual aminoacids in HEWL simulations.

References (1) Ninham, B. W.; Parsegian, V. A. Electrostatic potential between surfaces bearing ionizable groups in ionic equilibrium with physiologic saline solution. J. Theor. Biol. 1971, 31, 405–428. (2) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Am. Chem. Soc. 1990, 112, 6127–6129. (3) Warwicker, J.; Watson, H. C. Calculation of the electric potential in the active site cleft due to α-helix dipoles. J. Mol. Biol. 1982, 157, 671–679. (4) Bashford, D.; Karplus, M. pKa ’s of ionizable groups in proteins: atomic detail from a continuum electrostatic model. Biochemistry-US 1990, 29, 10219–10225. (5) Sharp, K. A.; Honig, B. Electrostatic Interaction in Macromolecules: theory and applications. Annu. Rev. Biophys. Bio. Chem. 1990, 19, 301–332. (6) Yang, A.-S.; Gunner, M. R.; Sampogna, R.; Sharp, K.; Honig, B. On the calculation of pKa s in proteins. Proteins: Struct. Funct. Genet. 1993, 15, 252–265. (7) Simonson, T. Macromolecular electrostatics: continuum models and their growing pains. Curr. Opin. Struct. Biol. 2001, 11, 243–252.

23



(8) Bashford, D. Macroscopic electrostatic models for protonation states in proteins. Front. Biosci 2004, 9, 1082–1099. (9) Rabenstein, B.; Knapp, E.-W. Calculated pH-dependent population and protonation of carbon-monoxy-myoglobin conformers. Biophys. J. 2001, 80, 1141–1150. (10) Wang, L.; Zhang, M.; Alexov, E. DelPhiPKa web server: predicting pKa a of proteins, RNAs and DNAs. Bioinformatics 2016, 32, 614–615. (11) Antosiewicz, J.; Porschke, D. The nature of protein dipole moments: experimental and calculated permanent dipole of alpha-chymotrypsin. Biochemistry-US 1989, 28, 10072–10078. (12) Beroza, P.; Fredkin, D. R.; Okamura, M. Y.; Feher, G. Protonation of interacting residues in a protein by a Monte Carlo method: application to lysozyme and the photosynthetic reaction center of Rhodobacter sphaeroides. Proc. Natl. Acad. Sci. USA 1991, 88, 5804–5808. (13) Baptista, A. M.; Soares, C. M. Some Theoretical and Computational Aspects of the Inclusion of Proton Isomerism in the Protonation Equilibrium of Proteins. J. Phys. Chem. B 2001, 105, 293–309. (14) Teixeira, V. H.; Cunha, C. C.; Machuqueiro, M.; Oliveira, A. S. F.; Victor, B. L.; Soares, C. M.; Baptista, A. M. On the Use of Different Dielectric Constants for Computing Individual and Pairwise Terms in Poisson-Boltzmann Studies of Protein Ionization Equilibrium. J. Phys. Chem. B 2005, 109, 14691–14706. (15) Mertz, J. E.; Pettitt, B. M. Molecular dynamics at a constant pH. Int. J. High Perform. C. 1994, 8, 47–53. (16) Baptista, A. M.; Martel, P. J.; Petersen, S. B. Simulation of protein conformational

24


Page 24 of 31

Page 25 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


freedom as a function of pH: constant-pH molecular dynamics using implicit titration. Proteins Struct. Funct. Bioinf. 1997, 27, 523–544. (17) Börjesson, U.; H¨ unenberger, P. H. Explicit-solvent molecular dynamics simulation at constant pH: methodology and application to small amines. J. Chem. Phys. 2001, 114, 9706. (18) Baptista, A. M.; Teixeira, V. H.; Soares, C. M. Constant-pH molecular dynamics using stochastic titration. J. Chem. Phys. 2002, 117, 4184–4200. (19) Burgi, 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. (20) Dlugosz, M.; Antosiewicz, J. M. Constant-pH molecular dynamics simulations: a test case of succinic acid. Chem. Phys. 2004, 302, 161–170. (21) Lee, M. S.; Salsbury, F. R.; Brooks III, C. L. Constant-pH molecular dynamics using continuous titration coordinates. Proteins Struct. Funct. Bioinf. 2004, 56, 738–752. (22) Mongan, J.; Case, D. A.; McCammon, J. A. Constant pH molecular dynamics in generalized Born implicit solvent. J. Comput. Chem. 2004, 25, 2038–2048. (23) Machuqueiro, M.; Baptista, A. M. Constant-pH Molecular Dynamics with Ionic Strength Effects: Protonation–Conformation Coupling in Decalysine. J. Phys. Chem. B 2006, 110, 2927–2933. (24) Stern, H. A. Molecular simulation with variable protonation states at constant pH. J. Chem. Phys. 2007, 126, 164112. (25) Donnini, S.; Tegeler, F.; Groenhof, G.; Grubm¨ uller, H. Constant pH molecular dynamics in explicit solvent with λ-dynamics. J. Chem. Theory Comput. 2011, 7, 1962–1978.

25



(26) Itoh, S. G.; Damjanović, A.; Brooks, B. R. pH replica-exchange method based on discrete protonation states. Proteins Struct. Funct. Bioinf. 2011, 79, 3420–3436. (27) Wallace, J. A.; Shen, J. K. Continuous constant pH molecular dynamics in explicit solvent with pH-based replica exchange. J. Am. Chem. Soc. 2011, 7, 2617–2629. (28) 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. (29) 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. (30) 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. (31) Bernardi, R. C.; Melo, M. C. R.; Schulten, K. Enhanced sampling techniques in molecular dynamics simulations of biological systems. BBA-Gen. Subjects 2015, 1850, 872–877. (32) Christen, M.; van Gunsteren, W. F. On searching in, sampling of, and dynamically moving through conformational space of biomolecular systems: A review. J. Comput. Chem. 2008, 29, 157–166. (33) Morrow, B. H.; Koenig, P. H.; Shen, J. K. Self-Assembly and Bilayer–Micelle Transition of Fatty Acids Studied by Replica-Exchange Constant pH Molecular Dynamics. Langmuir 2013, 29, 14823–14830. (34) Williams, S. L.; De Oliveira, C. A. F.; McCammon, J. A. Coupling constant pH molecular dynamics with accelerated molecular dynamics. J. Chem. Theory Comput. 2010, 6, 560–568. 26


Page 26 of 31

Page 27 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(35) Bennett, W. F. D.; Chen, A. W.; Donnini, S.; Groenhof, G.; Tieleman, D. P. Constant pH simulations with the coarse-grained MARTINI model - Application to oleic acid aggregates. Can. J. Chemistry 2013, 91, 839–846. (36) Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 1999, 314, 141–151. (37) Sugita, Y.; Kitao, A.; Okamoto, Y. Multidimensional replica-exchange method for free-energy calculations. J. Chem. Phys. 2000, 113, 6042–6051. (38) 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. (39) Ben-Naim, A. Cooperativity and regulation in biochemical processes; Kluwer Academic/Plenum Publishers, 2001. (40) Haynes, W. M.; Lide, D. R.; Bruno, T. J. Handbook of chemistry and physics; CRC Press, 2017. (41) Bartik, K.; Redfield, C.; Dobson, C. M. Measurement of the individual pKa values of acidic residues of hen and turkey lysozymes by two-dimensional 1H NMR. Biophys. J. 1994, 66, 1180–1184. (42) 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. (43) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of state calculations by fast computing machines. J. Chem. Phys. 1953, 21, 1087–1092. 27



(44) Hritz, J.; Oostenbrink, C. Hamiltonian replica exchange molecular dynamics using soft-core interactions. J. Chem. Phys. 2008, 128, 144121. (45) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press, USA, 1987; pp 191–198. (47) Vila-Vi¸cosa, D.; Silva, T. F. D.; Slaybaugh, G.; Reshetnyak, Y. K.; Andreev, O. A.; Machuqueiro, M. Membrane-Induced p K a Shifts in wt-pHLIP and Its L16H Variant. J. Chem. Theory Comput. 2018, 14, 3289–3297. (48) Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656–1676. (49) Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; Van Gunsteren, W. F. Definition and testing of the GROMOS force-field versions 54A7 and 54B7. Eur. Biophys. J. 2011, 40, 843–856. (50) Reis, P. B. P. S.; Vila-Vi¸cosa, D.; Campos, S. R. R.; Baptista, A. M.; Machuqueiro, M. Role of Counterions in Constant-pH Molecular Dynamics Simulations of PAMAM Dendrimers. ACS Omega 2018, 3, 2001–2009. (51) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Postma, J. P. M. A Consistent Empirical Potential for Water-Protein Interactions. Biopolymers 1984, 23, 1513–1518. (52) Tironi, I. G.; Sperb, R.; Smith, P. E.; van Gunsteren, W. F. A generalized reaction field method for molecular dynamics simulations. J. Chem. Phys. 1995, 102, 5451–5459.

28


Page 28 of 31

Page 29 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(53) Heinz, T. N.; van Gunsteren, W. F.; H¨ unenberger, P. H. Comparison of four methods to compute the dielectric permittivity of liquids from molecular dynamics simulations. J. Chem. Phys. 2001, 115, 1125–1136. (54) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (55) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182–7190. (56) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. Rapid grid-based construction of the molecular surface and the use of induced surface charge to calculate reaction field energies: Applications to the molecular systems and geometric objects. J. Comput. Chem. 2002, 23, 128–137. (57) Gilson, M. K.; Sharp, K. A.; Honig, B. Calculating the Eletrostatic Potential of Molecules in Solution: Method and Error Assessment. J. Comput. Chem. 1987, 9, 327–335. (58) Machuqueiro, M.; Baptista, A. M. Acidic range titration of HEWL using a constant-pH molecular dynamics method. Proteins Struct. Funct. Bioinf. 2008, 72, 289–298. (59) Machuqueiro, M.; Baptista, A. M. Is the prediction of pKa values by constant-pH molecular dynamics being hindered by inherited problems? Proteins Struct. Funct. Bioinf. 2011, 79, 3437–3447. (60) Aune, K. C.; Tanford, C. Thermodynamics of the denaturation of lysozyme by guanidine hydrochloride. I. Dependence on pH at 25◦ . Biochemistry-US 1969, 8, 4579–4585. (61) Aune, K. C.; Tanford, C. Thermodynamics of the denaturation of lysozyme by guanidine hydrochloride. II. Dependence on denaturant concentration at 25◦ . Biochemistry-US 1969, 8, 4586–4590.

29



(62) Hritz, J.; Oostenbrink, C. Optimization of replica exchange molecular dynamics by fast mimicking. J. Chem. Phys. 2007, 127, 204104. (63) Teixeira, V. H.; Vila-Vi¸cosa, D.; Reis, P. B. P. S.; Machuqueiro, M. pKa Values of Titrable Amino Acids at the Water/Membrane Interface. J. Chem. Theory Comput. 2016, 12, 930–934.

30


Page 30 of 31

Page 31 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical TOC Entry

31


A pH Replica Exchange scheme in the Stochastic Titration constant

Recommend Documents