Does Replica Exchange with Solute Tempering ... - ACS Publications

Article pubs.acs.org/JCTC

Does Replica Exchange with Solute Tempering Efficiently Sample Aβ Peptide Conformational Ensembles? Amy K. Smith, Christopher Lockhart, and Dmitri K. Klimov* School of Systems Biology and Computational Materials Science Center, George Mason University, Manassas, Virginia 20110, United States S Supporting Information *

ABSTRACT: We have applied replica exchange with solute tempering (REST) molecular dynamics to study a short fragment of the Aβ peptide, Aβ25-35, in water and a much larger system incorporating two Aβ10-40 peptides binding to the zwitterionic dimyristoylphosphatidylcholine (DMPC) bilayer. As a control, we used traditional replica exchange molecular dynamics (REMD) applied to the same systems. Our objective was to assess the practical utility of REST simulations. Taken together, our results suggest four conclusions. First, compared to REMD, the number of replicas in REST simulations can be reduced four to five times without affecting the temperature range or compromising an efficient random walk of REST replicas over temperatures. Second, although overall REST produces much fewer conformational states than REMD, there is no substantial difference in the collection of unique states for the wild-type replica in REST and REMD, especially for the system featuring Aβ peptides binding to the lipid bilayer. Third, we performed a thorough comparison of REST and REMD equilibrium conformational ensembles, including thermal averages and probability distributions. REST reproduces REMD data extremely well for the system of Aβ peptides binding to the DMPC lipid bilayer. The agreement between REST and REMD equilibrium sampling of Aβ25-35 in water is less perfect, but it improves with addition of new REST simulations. Surprisingly, REST demonstrates much better convergence for the system featuring ordered peptides binding to lipid bilayer rather than for a small unstructured peptide solvated in water. Fourth, REST delivers its full computational advantage over REMD when applied to peptides interacting with lipid bilayers. For peptides solvated in water, REST does not appear to offer computational gain but may make replica simulations practically feasible due to a lower requirement for parallel computing environments. Our study is expected to facilitate wider application of REST in biomolecular simulations.

1. INTRODUCTION Molecular dynamics (MD) simulations have become an important tool probing the mechanisms of biological processes at all-atom resolution. However, most biomolecular systems are characterized by rugged free-energy landscapes featuring deep minima surrounded by high barriers. Consequently, straightforward application of MD, in which external conditions and interaction parameters are kept fixed during simulation, typically results in a biased conformational sampling due to trapping of biomolecules in local free-energy minima. Lack of exhaustive conformational sampling makes it difficult to compare MD data with experimental observations reflecting bulk averages. In an attempt to overcome these issues, replica exchange molecular dynamics (REMD) simulations1−3 were introduced. Traditional (canonical or NVT) REMD simulations enhance conformational sampling through replication of a molecular system R times and assigning to each replica a unique temperature. In typical REMD implementations, the replica with the lowest temperature corresponds to normal (wild-type) conditions, whereas other R − 1 replicas are assigned elevated, exponentially distributed temperatures intended to accelerate conformational sampling.4 Periodically, replicas adjacent along the temperature scale are attempted to exchange. The © XXXX American Chemical Society

probabilities of accepting such exchanges are derived from the detailed balance condition, which preserves the Boltzmann probability distributions at individual temperatures. In practice, exchange attempts are governed by Metropolis Monte Carlo simulations in the replica space. If the REMD ensemble is simulated for a sufficiently long time, the replicas execute a random walk over temperatures, allowing them to escape local minima via simulations at elevated temperatures. An advantage of REMD is its flexibility in employing different conditions, which might enhance conformational sampling. For example, the REMD algorithm has been extended to isobaric−isothermal (NPT)5 or surface tension (NPγT)6 ensembles. Furthermore, a special version of REMD termed Hamiltonian is based on creating replicas with scaled atomic interaction parameters.7 Because scaling down these parameters promotes conformational transitions, a random walk of replicas over the range of interaction parameters prevents their trapping in local freeenergy minima. In the past, REMD has been widely applied to study the thermodynamics of protein folding, aggregation, and Received: July 2, 2016

A

DOI: 10.1021/acs.jctc.6b00660 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation binding.8−15 REMD has also been used extensively by our group16−18 to study equilibrium behavior of Aβ peptides. Practical use of REMD is restricted by severe computational requirements.4,19 To achieve equilibrium conformational sampling, one needs to consider a large number of replicas, R, covering a wide temperature range. Moreover, the probability of exchanging the replicas depends on the relative standard deviation of potential energy, which scales as N−1/2 with the number of atoms, N. Therefore, to ensure sufficient mixing of replicas across the temperature scale, i.e., to reach the average exchange probability of 0.2−0.3,4,20 it is necessary to use small intervals between adjacent temperatures, further increasing the number of required replicas. The need for a large number of replicas can be computationally restrictive, as expensive parallel computing environments are needed to conduct REMD simulations.19,21 To lower REMD computational costs, Liu et al.22 have developed a new version of Hamiltonian REMD termed replica exchange with solute tempering (REST) that scales the contribution of solute−solvent and solvent−solvent energies to the total energy function. In REST, the scaling factors depend on replica temperatures and are deliberately selected to strengthen solvent interactions at elevated temperatures. As a result, REST keeps the solvent effectively “cold” in all replicas, whereas the solute is simulated at different temperatures as in traditional REMD. Note that REST always keeps at least one energy function intact or wild-type, whereas others are affected by scaling. Such targeted scaling of molecular interactions involving solvent eliminates their contribution to exchange probabilities. Consequently, the replica exchange probabilities depend exclusively on the contribution from solute atoms, which constitute a small fraction of all atoms. Because distributions of solute energies are relatively broad, few replicas are needed to cover the desired temperature range and achieve sufficiently high probabilities of exchanges. Published studies and our own estimates have shown that REST can reduce the number of replicas 5- to 10-fold, thus considerably alleviating REMD computational requirements. Use of REST has gained traction in recent years,23,24 especially in ligand binding,25−27 and it has been implemented in several popular MD programs.28,29 However, despite potential advantages over traditional REMD, the REST efficiency has not been systematically studied, nor has it been applied to complex biomolecular systems involving, for example, binding of peptides to lipid bilayers. These systems pose a special challenge for REST compared to lipid-free aqueous simulations because, to take advantage of a small number of replicas, REST needs to include lipids in a “solvent” category in addition to water. However, cold lipids may hamper peptide conformational sampling, thus rendering REST inefficient. It is also conceivable that, since the overall scope of conformational sampling produced by REST is limited compared to REMD, the diversity of conformations collected with the wild-type REST energy does not match the scope of sampling produced by REMD for the low temperatures. As a result, much longer REST simulations might be needed to reach parity with REMD. To explore the issues posed above, we compared the performance of two sampling algorithms, REMD and REST, by applying them to two biomolecular systems. First, we used REMD and REST to study the behavior of a small Aβ fragment, Aβ25-35, in water. This fragment forms a random coil/turn conformation in a 20:80 hexafluoroisopropanol:water mixture30 and is amyloidogenic.31 Second, we considered a much larger

system incorporating two Aβ10-40 peptides binding to the zwitterionic dimyristoylphosphatidylcholine (DMPC) bilayer. REMD simulations of this system have been reported by us earlier.32,33 By performing REST simulations for the same system, we compared the two sampling methods. In general, our goals in comparing REMD against REST are 2-fold. First, we check the ability of REST to provide exhaustive unbiased conformational sampling comparable with REMD. Second, we evaluate practical computational advantages offered by REST.

2. MODEL AND SIMULATION METHODS 2.1. Replica Exchange with Solute Tempering. Traditional REMD replicates a molecular system R times and assigns unique temperatures to each replica. Periodically, structures at neighboring temperatures are attempted to exchange following the detailed balanced condition, which preserves the Boltzmann distribution at each temperature. REST22 follows the same basic scheme, although its replicas differ not only with respect to their temperatures but also with respect to scaling factors applied to solute−solvent and solvent−solvent energies. In this sense, REST simulations are akin to Hamiltonian REMD,7 where each replica samples a different Hamiltonian. Following the scaling factors proposed by Camilloni et al.34 and later used by Best et al.,35 the enthalpy Hr of the isobaric−isothermal (NPT) system with the pressure P and volume V is defined as Hr = Ep +

β0 βr

Epw +

β0 βr

Ew + PV (1)

Here, Ep, Epw, and Ew are the solute−solute, solute−solvent, and solvent−solvent potential energies and the factor βr = (kBTr)−1 corresponds to the temperatures with the indices 0 ≤ r ≤ R − 1. In total, there are R different temperatures. Importantly, at the lowest temperature r = 0, the energy function is unperturbed or wild-type. All other replicas simulated at the temperatures r > 0 feature perturbed Hamiltonians and serve to enhance conformational sampling. In REST, force constants for solvent atom bond lengths, angles, and dihedrals as well as Lennard-Jones well depths, ϵ, are multiplied by the factors β0/βr. Solvent partial charges are multiplied by the factors β0 /βr . These scaling factors are applied within the MD source code to minimize loss of numerical accuracy as it has been implemented in GROMACS28 and NAMD.29 REST simulations attempt to swap the neighboring replicas r and r + 1 with the system configurations Xr and Xr+1, respectively. Then, acceptance of a replica exchange follows from the Metropolis criterion with the probability ω = min[1, e−Δ], where Δ is Δ = βr (Hr(X r + 1) − Hr(X r )) + βr + 1(Hr + 1(X r ) − Hr + 1(X r + 1)) (2)

Taking into account the scaling of solvent−solvent terms in eq 1, their contribution to the exchange probability, ω, is eliminated and Δ is reduced to Δ = (βr − βr + 1)(Ep(X r + 1) − Ep(X r )) + ( β0βr −

β0βr + 1 )(Epw (X r + 1) − Epw (X r ))

+ (βr − βr + 1)P(Vr + 1 − Vr ) B

(3)


Article

Journal of Chemical Theory and Computation

Simulations were performed using the MD program NAMD38 updated with REST support.29 All simulations utilized periodic boundary conditions. Electrostatic interactions were computed using Ewald summation, whereas van der Waals interactions were smoothly switched off in the interval from 8 to 12 Å. Covalent bonds associated with hydrogen atoms were constrained by the SHAKE algorithm. Temperature was controlled with underdamped Langevin simulations of “virtual” solvent with the damping coefficient γ = 5 ps−1. Pressure was held constant at 1 atm by applying the Langevin piston method. The integration step was set to 1 fs. Replica exchanges were attempted every 2 ps between neighboring replicas along the temperature scale. The first, small system shown in Figure 1a included a single Aβ25-35 peptide in a water box with the initial dimensions 36.7 Å × 36.7 Å × 36.7 Å containing 1602 water molecules. REMD simulations serving as control featured R = 16 replicas with the temperatures distributed exponentially from 320 to 400 K. REST simulations required 4-fold fewer replicas (R = 4) with the temperatures distributed exponentially in the same range. Because in REST simulations solute and solvent must have a zero net charge, a single chloride counterion was included as solute along with the Aβ25-35 peptide. Five REMD and nine REST independent trajectories were produced. Each replica in a trajectory was simulated for 20 ns. In total, REMD accumulated 1.6 μs of sampling, whereas 0.72 μs was collected for REST. Initial unequilibrated portions of the trajectories were discarded, resulting in 1.2 and 0.63 μs of equilibrated sampling, respectively. The second, large simulation system displayed in Figure 1b contained two Aβ10-40 peptides co-incubated with 98 DMPC lipids (49 lipids per leaflet). The simulation box had the initial dimensions 55.6 Å × 55.6 Å × 81.4 Å containing 4356 water molecules. As a control, we used our previous REMD simulations of this system,32,33 which involved R = 40 replicas distributed exponentially from 320 to 430 K. In contrast to REMD, REST utilized 5-fold fewer replicas (R = 8) distributed exponentially from 330 to 430 K. REST simulations were performed using the greedy replica exchange algorithm,19 which improves the efficiency of REMD simulations by removing strict synchronization across all R replicas at each iteration. A semi-isotropic coupling scheme was used that couples the x and y box dimensions while allowing the z dimension to fluctuate independently. Aβ10-40 peptides and two sodium counterions were considered as a solute, whereas water and lipids were treated as solvent. Inclusion of lipids as solvent during REST simulations keeps them effectively cold, preventing bilayer disintegration at high REST temperatures. However, to facilitate comparison between REMD and REST, two constraints as implemented in our previous REMD simulations32 were also applied to REST simulations. Their purpose was to prevent bilayer disintegration at high REMD temperatures and peptide aggregation (see the Supporting Information). For both REMD and REST simulations, five independent trajectories were produced. In the initial structures, four or five amino acids (≤16% of all) were bound to the DMPC surface, but none were inserted; i.e., they did not occur below phosphorus atoms. Each replica was simulated for 20 ns. In all, REMD accumulated 4 μs of data, whereas REST collected 0.8 μs of sampling. Initial portions of the trajectories were discarded bringing the equilibrated sampling time to 3.32 and 0.66 μs, respectively. Due to the inclusion of two non-

From a practical viewpoint, it is easier to use the exchange criterion based on the total system enthalpy specified by eq 2 rather than eq 3. Additional details concerning our REST implementation are available in the Supporting Information. 2.2. Replica Exchange Simulations. Isobaric−isothermal (NPT) REMD and REST simulations of the Aβ25-35 peptide in water and of Aβ10-40 peptides bound to the DMPC bilayer were performed (Figure 1). Aβ25-35 simulations employed a

Figure 1. (a) Small simulation system consisting of Aβ25-35 peptide in water. (b) Large simulation system containing two Aβ10-40 peptides binding to the opposite leaflets of the DMPC bilayer. In both panels coloring reflects “hot” (peptide(s) in red) and “cold” (water and lipids in blue) REST partitions. Phosphorus atoms in panel b shown as large spheres fluctuate around ±zP. The bilayer midplane corresponds to z = 0 Å.

newer version of the CHARMM force field, CHARMM36.35 To preserve consistency with our previous study32 serving as a control, REST simulations of Aβ10-40 peptides bound to the DMPC bilayer utilized the CHARMM22 force field with CMAP corrections36 to model peptides and CHARMM3637 to model lipids. Both simulation systems used the CHARMMmodified TIP3P water model. Aβ peptides were simulated at neutral pH with acetylated and aminated termini. C


Article


than 6.5 Å. Aβ10-40 peptides bound to the DMPC bilayer were characterized by the probability distribution P(z;i), which maps the location of each amino acid i along the z axis normal to the bilayer (Figure 1b). Peptide−lipid contacts were computed by dividing the lipid molecule into five segments, choline (G1), phosphate (G2), glycerol (G3), and two fatty acid tails (G4 and G5) and applying the same definition as for an intrapeptide side chain contact. The number density nl(r,z) of heavy lipid atoms was computed as a function of the distance r to the peptide center of mass in the (x, y) plane and the distance z to the bilayer midplane (Figure 1b). Last, we computed the carbon− deuterium order parameter, SCD, to probe ordering of lipid tails defined as

interacting Aβ peptides, the effective amount of peptide sampling is doubled. 2.3. Measures of Sampling. To compare REMD and REST algorithms, we investigated their performance and quality of conformational sampling using four measures. First, we computed the enthalpy H probability distributions, P(H;T), at each simulation temperature, T. Sufficient overlap of adjacent P(H;T) is a prerequisite for efficient replica exchange.4,20 Second, we computed the replica mixing parameter, m(T), introduced by Han and Hansmann,39 R

m (T ) = 1 −

∑r = 1 tr 2 R

∑r = 1 tr

(4)

where tr denotes the total number of iterations replica r spends at temperature T. The parameter m(T) was computed for each temperature T utilized in REST and REMD. If all R replicas are equally represented across the entire temperature distribution, then m(T) = 1 − 1/√R, which is a temperature independent theoretical optimum value of m(T). Approach of m(T) to its optimal value implies a random walk of replicas over temperatures. Third, we computed the number of unique states, Ns(τsim), sampled during simulations, where τsim is the equilibrium simulation time accumulated by all replicas r or solely at the lowest temperature (r = 0). A unique state (H,X) is defined by an enthalpy H and structural quantity X (see Results and the Supporting Information). The bin size for H was set to 2 kcal/ mol. The purpose of computing Ns(τsim) is to investigate the ability of REMD or REST to exhaust conformational space. Saturation of Ns(τsim) constitutes a necessary but not sufficient condition for simulation convergence. Alternatively, following other studies40 we could cluster peptide conformations and count unique clusters as a function of τsim. However, in our simulations this approach would not directly account for peptide−lipid interactions. Fourth, we considered the sampling errors in REMD and REST. The difference between REST and REMD sampling was assessed by computing the relative rootmean-square deviation (rRMSD) defined for the residuespecific quantity X(i) as

SCD =

3. RESULTS To evaluate the utility of REST simulations, we considered a small system consisting of the Aβ25-35 peptide in water (4,968 atoms) and a much larger system, which includes two Aβ10-40 peptides binding to the DMPC bilayer (25,586 atoms). For both systems we have performed REST simulations and used as a control standard REMD simulations. Our objectives are (i) to compare the quality and efficiency of REST and REMD conformational sampling and (ii) to test the consistency between the equilibrium conformational ensembles produced by REMD and REST. 3.1. Comparing REST and REMD Sampling. Aβ25-35 Peptide in Water. We first considered the enthalpy probability distributions, P(βH;T), at each of the four REST temperatures. Figure 2a confirms that these distributions show significant overlaps resulting in an average exchange rate of 31%. The distributions P(H;T) obtained from the control REMD simulations are presented in Figure S1a for all 16 replicas characterized by an average exchange rate of 29%. Importantly, the REMD probability distributions of enthalpies collected at the temperatures of REST simulations are clearly separated featuring virtually no overlaps. To evaluate the random walk of replicas over temperatures, we have computed the replica mixing parameter m(T) at each

((1/N ) ∑i ( XREST(i) − XREMD(i) )2 )1/2 N

(1/N ) ∑i

(6)

where θ is the angle between the bilayer normal and a C−H bond in the sn-2 fatty acid tail. SCD was computed twice, once for lipids proximal to Aβ and a second time for lipids distant from Aβ. Following our previous study,33 proximal lipids were defined as those occurring within the distance Rc from the center of mass of Aβ along the (x, y) plane. Distant lipids were those outside this radius. The specific value of Rc corresponds to the distance r where peptide−lipid contacts reach maximum. For REST simulations, Rc = 13 Å, whereas for REMD Rc = 14 Å.33 Thermodynamic averages of structural quantities were computed using the weighted histogram analysis method (WHAM)42−44 as described in the Supporting Information. WHAM was used to compute the thermodynamic averages at the lowest temperature with a wild-type Hamiltonian. For the Aβ25-35 fragment in water, this temperature was 320 K, whereas for Aβ10-40 peptides bound to the DMPC bilayer it was 330 K. A slightly higher temperature for Aβ10-40 was chosen to match our previous simulations32 and to ensure that the DMPC bilayer samples a liquid disordered phase. Unless noted otherwise, all thermodynamic quantities are reported at these temperatures using all available trajectories.

N

rRMSD(X ) =

3 cos2 θ − 1 2

XREMD(i) (5)

where N is the number of amino acids. rRMSD evaluates the difference in REST and REMD thermodynamic averages denoted as ⟨...⟩ with respect to the average amplitude of X. Additional comparison of REST and REMD is offered by the error metric, ρ, defined as a ratio of the root-mean-square deviation (RMSD) between REST and REMD thermodynamic averages of a quantity X and the REMD sampling error δX. If X is residue-specific, then ρ = RMSD(X)/δX, where RMSD(X) = ((∑Ni=1(⟨XREST(i)⟩ − ⟨XREMD(i)⟩)2)/N)1/2 and the average sampling error δX = (∑Ni=1 δX(i))/N. To characterize peptides and lipids, we used the following structural quantities. Aβ secondary structure was computed using the program STRIDE.41 Amino acids were considered in a helical state if they sampled α-, 310-, or π-helices. We have also computed the probability distributions P(ϕ,ψ) of peptide backbone ϕ and ψ dihedral angles. Intrapeptide contacts between two amino acids i and j were considered formed when the distance between their heavy atom centers of mass was less D


Article


Figure 2. Probability distributions of enthalpies P(βH;T) in REST simulations for Aβ25-35 (a) and Aβ10-40 (b) systems. Colors distinguish four or eight REST temperatures used for Aβ25-35 and Aβ10-40. Significant overlap between adjacent P(βH;T) distributions leads to frequent replica exchanges between temperatures. Following the REST formalism, enthalpies H are multipled by the Boltzmann factor β.

Figure 3. Replica mixing parameter m(T) computed as a function of REMD and REST temperatures T for Aβ25-35 (a) and Aβ10-40 (b) systems. REST and REMD data are represented by black and red circles. The dashed lines mark the optimal theoretical values of m(T) for REST simulations. REMD values of m(T) are adjusted by subtracting 1/ RREST − 1/ RREMD , where RREST and RREMD are the numbers of REST and REMD replicas for a given system. Close agreement between computed and optimal values of m(T) indicates efficient mixing of replicas across temperatures.

temperature, T, simulated in REMD and REST (see Model and Simulation Methods). Figure 3a presents m(T) averaged over all REMD and REST trajectories. According to eq 4, 16 REMD and four REST replicas would result in optimal theoretical values of 0.75 and 0.50, respectively. In the figure, near-constant values of m(T) almost coincide with their optimal values for both simulations indicating near-ideal mixing of replicas over temperatures. The distribution of replicas over temperatures during REST simulations is directly visualized in Figure S2a. Consistent with the behavior of m(T), Figure S2a reveals a random walk of replicas with no indications of their trapping at any temperature. The convergence of REMD and REST sampling can be tested by computing Ns of unique states (H, X) sampled at least once during equilibrium simulation time (see Model and Simulation Methods). As a structural quantity X, we chose the total number of intrapeptide contacts, Cp. Figure 4a shows Ns as a function of equilibrium simulation time, τsim, collected at all REST or REMD temperatures. It is seen that the increase in Ns

dramatically slows if τsim ≳ 0.1 μs for REST or ≳0.4 μs for REMD. This behavior of Ns signals a gradual exhaustion of new states that is a prerequisite for sampling convergence. More strikingly, the figure underscores a vast gap in the number of unique states sampled by REST and REMD, which is a consequence of a 4-fold difference in the numbers of replicas and different conditions employed in REST and REMD. The inset to Figure 4a compares Ns collected exclusively at the lowest temperature (320 K) with the wild-type energy function. It again shows that REST gradually exhausts new states at this temperature. Another observation is that REST visits a slightly smaller (∼10%) number of unique states than REMD at large τsim. In fact, the REST equilibrium simulation time at the lowest temperature must be roughly doubled to reach parity in the number of unique states with REMD. As a further measure of simulation convergence, we computed the sampling errors as a function of the number of E


Article


⟨T(i)⟩, for amino acids i (see Model and Simulation Methods) was computed as a function of the number of REST or REMD trajectories, M. Similar to the error in side chain interactions, Figure 5b reveals a decrease in δT with M for both algorithms. However, the figure also exhibits strikingly larger errors observed in REST compared to REMD. For example, for M = 5 the REMD error δT is ∼0.02, but the corresponding REST error is almost twice larger. It is seen that even when nine REST trajectories are produced, the REST turn error is still larger than the REMD error at M = 5. Simple extrapolation of δT(M) for REST suggests that it becomes equal to REMD δT(M=5) when M ≳ 20. We have also considered the error in peptide end-to-end distance, r1N. However, we found that the REST sampling errors in r1N are lower than in REMD for a given number of trajectories M. Thus, our analysis implies that for some structural quantities (such as turn propensities) the number of REST trajectories must be considerably increased (e.g., 4-fold) compared to REMD to reach equal errors. Aβ10-40 Peptides Binding to DMPC Bilayer. To further investigate the efficiency and quality of the REST algorithm, we studied a large system consisting of two Aβ10-40 peptides binding to the DMPC bilayer. As for Aβ25-35, we first compared the enthalpy distributions computed from REST and REMD simulations. Figure 2b shows the enthalpy probability distributions P(βH;T) for the eight REST temperatures, whereas Figure S1b presents the distributions P(H;T) for the 40 REMD temperatures. Both figures demonstrate significant overlap between adjacent enthalpy distributions, which is a prerequisite for efficient replica exchange. Consequently, the average replica exchange rates were 22% and 24% for REST and REMD, respectively. Figure S1b also indicates that the REMD enthalpy distributions computed at eight REST temperatures, which show significant overlap in REST, have virtually no overlap in REMD and thus would not perform replica exchanges. REST simulations were also evaluated by computing the replica mixing parameter, m(T). According to Model and Simulation Methods, when the number of REST replicas R is 8, the optimum value of m(T) is 0.65, whereas increasing in REMD R to 40 raises the optimum value to 0.84. Figure 3b compares m(T) computed for REST and REMD simulations with their optimum values and demonstrates a good agreement between optimum and calculated values for most temperatures. Indeed, the average calculated values of m(T) for REST and REMD are 0.63 and 0.82, which are close to the respective optimum values. Such agreement indicates an efficient mixing of replicas resulting in their approximately equal distribution over temperatures. A similar conclusion can be drawn from Figure S2b, which depicts the random walk of replicas over temperatures for a typical REST trajectory. Next, we describe the number of unique states sampled during REST or REMD simulations. A state (H,X) is specified by the enthalpy H and a structural quantity X, where X represents either Cp or Cl. Figure 4b shows the Ns for these two choices of X as a function of τsim collected at all REST and REMD temperatures. Both Ns computed using REMD or REST level off with simulation progress (for τsim ≳ 2.5 μs or ≳ 0.4 μs, respectively) indicating an exhaustion of new states. However, similar to Aβ25-35 simulations, a large gap between the final numbers of unique states collected in REMD and REST exists reflecting a dramatically larger scope of REMD simulations. It is important that the inset to Figure 4b demonstrates that REMD and REST sampling become nearly

Figure 4. (a) Number Ns of unique states (H, Cp) defined using the wild-type enthalpy H and the number of intrapeptide contacts Cp, which are sampled at least once during equilibrium simulation time τsim collected in all replicas. The data are computed for the Aβ25-35 system. (b) Ns (H, X) defined using the wild-type enthalpy H and structral quantity X collected in all replicas. The data are computed for the Aβ10-40 system. Solid and dashed lines represent Ns when X is given by the number of intrapeptide contacts Cp or the number of peptide−lipid contacts Cl. In both panels, REST and REMD data are shown in black and red. The insets present the Ns collected exclusively at the temperature of interest (320 K for Aβ25-35 and 330 K for Aβ10-40). Slowdown in Ns growth with simulation progress signifies exhaustion of new states.

REST or REMD trajectories, M, treating each trajectory as an independent sample. We first considered the average error δCp in the probabilities of Aβ25-35 intrapeptide contacts, ⟨Cp(i,j)⟩, between amino acids i and j (defined similarly to the error in residue-specific quantity X in Model and Simulation Methods). Figure 5a, which shows δCp as a function of M, indicates that the REST and REMD sampling errors monotonically decrease with M as expected. The figure also demonstrates that the REST errors are consistently larger than those of REMD, implying that to reach REMD error magnitude, one needs to produce at least one extra REST trajectory. To get further insight, we analyzed the average sampling errors in amino acid turn propensities. (Turn propensity was selected, because it represents the only ordered structure in Aβ25-35.) Specifically, the average error δT in the equilibrium turn propensities, F


Article


Figure 5. (a) Average sampling error δCp in the probabilities of side chain contacts in the Aβ25-35 peptide as a function of the number of REMD or REST trajectories M. (b) Average sampling error in residue-specific turn propensities δT as a function of the number of REMD or REST trajectories M for the Aβ25-35 peptide. (c) Average sampling error in residue-specific helix propensities δH as a function of the number of REMD or REST trajectories M for the Aβ10-40 peptide. (d) Average sampling error δz in the positions of Aβ10-40 amino acids along the bilayer normal as a function of the number of REMD and REST trajectories M. In all panels data for REMD and REST are in red and black. The figure reveals minor differences in REST and REMD errors for a bilayer-bound peptide but points out that for a small unstructured peptide in water REST generates larger errors than REMD.

identical when the collection of new states is restricted to the (wild-type) temperature of 330 K. Thus, for the Aβ10-40 system the REST algorithm collects the conformational states at 330 K in a way remarkably similar to REMD (see Discussion). Finally, we investigated the sampling errors of REST and REMD simulations. We first considered the average sampling error in the amino acid helix propensities, δH. This quantity was chosen because, according to our previous REMD simulations, the structural ensemble of Aβ10-40 bound to the zwitterionic bilayer is predominately helical.32 Figure 5c compares the helical sampling errors δH for REST and REMD. As M increases, the REST and REMD sampling errors remain close to each other and concomitantly decrease. We also plot in Figure 5d the error δz in the position ⟨z(i)⟩ of amino acid i along the bilayer normal averaged over all amino acids as a function of M. The quantity ⟨z(i)⟩ reflects the interactions between Aβ amino acids and lipids, which are kept cold in all REST replicas. Figure 5d demonstrates that, despite the inclusion of lipids in the cold category, the sampling errors δz in REST and REMD are nearly identical. 3.2. Comparing REST and REMD Conformational Ensembles. Aβ25-35 Peptide in Water. The second objective of our study is to test the consistency between equilibrium conformational ensembles produced by REMD and REST. As above, we first consider the Aβ25-35 peptide in water.

Computation of secondary structure sampled by REST and REMD shows that, consistent with the experiments,30,45 Aβ2535 adopts mostly random coil conformations augmented by significant turn propensity. Indeed, the average fractions of random coil structure, ⟨RC⟩, in REST and REMD simulations are identical (0.69 ± 0.03) indicating that in both conformational ensembles the average numbers of random coil amino acids are 7.6 ± 0.3. The REST and REMD fractions of turn states are 0.29 ± 0.03 and 0.26 ± 0.02 corresponding to the numbers of turn amino acids of 3.2 ± 0.3 and 2.9 ± 0.2. Occurrence of helix or β-sheet is generally rare. Thus, the REST- and REMD-generated secondary structure propensities differ by no more than about 10% and are equal within their sampling errors. Figure 6a takes a next step by comparing the equilibrium turn fractions, ⟨T(i)⟩, for individual amino acids i. Good agreement between REST and REMD is seen within the C-terminal (positions 30−35), while the N-terminal (25−29) reveals elevated turn fraction for REST. The relative RMSD between REST and REMD values of ⟨T(i)⟩ (see Model and Simulation Methods) is 0.20, reflecting generally small but not negligible differences between REST and REMD data. REMD and REST sampling is further compared by computing the error metric, ρ, which is the ratio of the RMSD between the REST and REMD thermodynamic averages of a quantity X and the REMD sampling error δX. Using as X the residue-specific ⟨T(i)⟩, we found that ρ = 2.7 suggesting that the difference G


Article


Aβ25-35 tertiary structure was probed by computing intrapeptide contacts. The contact map, ⟨Cp(i,j)⟩, which reports the equilibrium probabilities of forming contacts between amino acids i and j in REST simulations, is shown in Figure 7a.

Figure 7. (a) Contact map (upper half), ⟨Cp(i,j)⟩, reporting interactions between Aβ25-35 amino acids i and j in REST simulations. Difference contact map (lower half) ⟨ΔCp(i,j)⟩, which visualizes the changes in contact probabilities between REST and REMD. (b) Probability distributions, P(Rg), of radius of gyration Rg for Aβ25-35 peptide. Black and red bars represent REST and REMD data. The panels show that REST and REMD distributions are generally consistent.

Figure 6. (a) Equilibrium turn fractions ⟨T(i)⟩ for Aβ25-35 amino acids i observed in REST (in black) and REMD (in red) simulations. Sampling errors are shown by vertical bars. The horizontal dashed line marks ⟨T⟩ = 0.5. (b) Difference ΔP(ϕ,ψ) in the probability distributions for Aβ25-35 backbone dihedral angles ϕ and ψ produced by REST and REMD. Both panels demonstrate a good but not perfect agreement between REST and REMD data.

Consistent with the prevalence of random coil conformations, the contact map reveals a lack of stable (>0.5) intrapeptide interactions. Indeed, the five most frequent contacts listed in Table 1 are all local (|i − j| < 5) and at best marginally stable (≲0.3). All long-range interactions (|i − j| ≥ 5) are completely disrupted occurring with the probabilities not exceeding 0.1. Accordingly, the average total number of side chain contacts in the Aβ25-35 fragment ⟨Cp⟩ is just 3.9. More importantly, the REST contact map is consistent with the control REMD ⟨Cp(i,j)⟩. From REMD simulations, the average number of

between REMD and REST averages is larger than the REMD sampling errors. This result is consistent with a large sampling error in REST seen in Figure 5b and a difference in turn propensities in the N-terminal. This issue is further assessed in the Discussion. To provide detailed analysis of Aβ25-35 backbone conformations, we have computed the probability distributions, P(ϕ,ψ), of backbone dihedral angles ϕ and ψ produced by REST and REMD. Figure 6b compares the REST and REMD P(ϕ,ψ) by displaying the difference ΔP(ϕ,ψ) = PREST(ϕ,ψ) − PREMD(ϕ,ψ). In general, ΔP(ϕ,ψ) ≈ 0 with the exception of the α-helix region, where |ΔP(ϕ,ψ)| ≲ 0.02. For a quantitative comparison, we computed the average relative RMSD between REST and REMD P(ϕ,ψ) distributions, which is defined similarly to turn fraction above. It follows that rRMSD is fairly large (0.37 when using nonzero values of P(ϕ,ψ)), whereas ρ is 2.2. This finding together with the turn structure analysis suggest that REST- and REMD-generated secondary structure propensities are reasonably consistent, although not identical.

Table 1. Five Most Frequent Intrapeptide Contacts in Aβ2535 Peptidea

a

H

rank

REST

REMD

1 2 3 4 5

Ser26-Lys28 (0.30 ± 0.03) Gly33-Met35 (0.29 ± 0.02) Gly29-Ile31 (0.22 ± 0.02) Ala30-Ile32 (0.20 ± 0.02) Gly25-Lys28 (0.19 ± 0.04)

Ser26-Lys28 (0.32 ± 0.02) Gly33-Met35 (0.30 ± 0.05) Gly25-Asp27 (0.25 ± 0.01) Ile31-Gly33 (0.22 ± 0.03) Ala30-Ile32 (0.20 ± 0.02)

The probabilities of forming a contact are in parentheses. DOI: 10.1021/acs.jctc.6b00660 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation Table 2. Secondary Structure Propensities in Aβ10-40 Sequence Regionsa structure ⟨H⟩ ⟨T⟩ ⟨RC⟩ ⟨S⟩ a

S1 0.18 (0.14 0.45 (0.46 0.37 (0.40 0.00 (0.01

± ± ± ± ± ± ± ±

S2 0.07 0.08) 0.05 0.06) 0.03 0.04) 0.00 0.00)

0.11 (0.10 0.41 (0.39 0.47 (0.50 0.01 (0.01

± ± ± ± ± ± ± ±

S3 0.03 0.03) 0.03 0.03) 0.04 0.01) 0.01 0.01)

0.37 (0.39 0.44 (0.41 0.19 (0.20 0.00 (0.00

± ± ± ± ± ± ± ±

S4 0.04 0.05) 0.06 0.06) 0.02 0.02) 0.00 0.00)

0.61 (0.65 0.11 (0.10 0.27 (0.26 0.01 (0.01

± ± ± ± ± ± ± ±

0.05 0.05) 0.05 0.05) 0.01 0.01) 0.01 0.00)

Secondary structure propensities for REMD32 are in parentheses.

contacts, ⟨Cp⟩, is found to be 4.0; i.e., it deviates from the REST value by about 3%, and similar to REST, no stable contacts are observed. The five most frequent REMD contacts are given in Table 1, which demonstrates that REST and REMD share three most frequent contacts with equal (within errors) probabilities. If the 10 most frequent contacts are considered, then eight of them are shared between REMD and REST. As in REST, longrange contacts in REMD are completely disrupted. To provide a direct comparison, we include in Figure 7a the difference contact map ⟨ΔCp(i,j)⟩ = ⟨Cp(i,j)⟩REST − ⟨Cp(i,j)⟩REMD, which reports the variations in contact probabilities between REST and REMD. The absolute value of ⟨ΔCp(i,j)⟩ does not exceed 0.1, and the relative RMSD is 0.34. If the computations are restricted to the contacts with significant probabilities (⟨Cp(i,j)⟩ > 0.1 max(⟨Cp(i,j)⟩)), rRMSD hardly changes (0.30). ρ comparing REST and REMD samplings of amino acid contacts was found to be 1.7. Finally, Aβ25-35 structure was evaluated by computing the probability distribution of radius of gyration, P(Rg). The REST conformational ensemble has ⟨Rg⟩ = 11.4 ± 0.1 Å, which is in good agreement with REMD-produced ⟨Rg⟩ = 11.6 ± 0.2 Å. The REST and REMD probability distributions, P(Rg), shown in Figure 7b are also consistent. This conclusion is confirmed by the low rRMSD between the two distributions (0.21), whereas the error metric for the probability distributions P(Rg) is ρ = 1.7. The analysis of Cl− counterion binding to Aβ can be found in the Supporting Information. Aβ10-40 Peptides Binding to DMPC Bilayer. We next compare the REST- and REMD-produced equilibrium conformational ensembles of Aβ10-40 peptides binding to the DMPC bilayer. We first consider the Aβ secondary structure. Overall, the fractions of helix ⟨H⟩, turn ⟨T⟩, random coil ⟨RC⟩, and β-structure ⟨S⟩ observed in REST are 0.38 ± 0.02, 0.31 ± 0.02, 0.31 ± 0.01, and 0.01 ± 0.00, respectively. The corresponding quantities for REMD32 are 0.39 ± 0.03, 0.30 ± 0.03, 0.31 ± 0.01, and 0.01 ± 0.00, showing an excellent (within less than 3% difference) agreement between the two algorithms. In Table 2, we further analyze the secondary structure propensities by dividing the Aβ10-40 into four sequence regions, s: S1 (residues 10−16), S2 (residues 17−21), S3 (residues 22−28), and S4 (residues 29−40). Importantly, for all regions the REST and REMD secondary structure propensities are within each other’s error. The REST data thus confirm the findings in our previous work32 that the helical state dominates Aβ conformations bound to the zwitterionic DMPC bilayer. To compare residue-specific helical propensities produced by REST and REMD, we plot in Figure 8a the helix fractions, ⟨H(i)⟩, for amino acids i. The relative RMSD computed between the REST and REMD distributions ⟨H(i)⟩ is 0.10, indicating minor deviations between the two algorithms.

Figure 8. (a) Comparison of the helical propensities ⟨H(i)⟩ for Aβ1040 amino acids i produced by REST (in black) and REMD32 (in red). The plot illustrates that the REST and REMD data are within their sampling errors shown by vertical bars. The horizontal dashed line marks ⟨H⟩ = 0.5. (b) Contact map (upper half), ⟨Cp(i,j)⟩, displaying the probabilities of contacts between Aβ10-40 amino acids i and j in REST simulations; difference contact map (lower half), ⟨ΔCp(i,j)⟩, comparing REST and REMD32. The difference contact map illustrates an excellent agreement between REST and REMD.

ρ, which measures the variations between REST and REMD sampling with respect to REMD error, is 0.57, implying that the difference between the algorithms is actually smaller than the REMD sampling error (see Discussion). Overall, these data demonstrate a remarkable consistency between the secondary structure propensities generated by REST and REMD. Following the analysis of Aβ25-35, the tertiary fold of the Aβ10-40 peptide was examined by computing the intrapeptide contact maps, ⟨Cp(i,j)⟩ (Figure 8b). The average number of contacts, ⟨Cp⟩ formed in REST simulations is 28.5 ± 1.0, which is equal to the REMD value within the sampling error (28.4 ± 0.7). Restricting the computations to long-range contacts (|j − i| ≥ 5), we find the number of such contacts to be ⟨CLR p ⟩ = 11.5 I


Article

Journal of Chemical Theory and Computation ± 1.0 and 11.1 ± 0.6 for REST and REMD, respectively, which are again equal within their errors. The top five stable contacts for REST and REMD simulations are listed in Table 3. It Table 3. Five Most Frequent Intrapeptide Contacts in Aβ1040 Peptidea

a

rank

REST

REMD32

1 2 3 4 5

Gly33-Val36 (0.87 ± 0.09) Gly33-Gly37 (0.82 ± 0.09) Ser26-Lys28 (0.77 ± 0.07) Asp23-Lys28 (0.71 ± 0.08) Ala21-Val24 (0.63 ± 0.03)

Gly33-Val36 (0.88 ± 0.08) Gly33-Gly37 (0.83 ± 0.09) Asp23-Lys28 (0.79 ± 0.06) Ser26-Lys28 (0.76 ± 0.08) Ala21-Val24 (0.74 ± 0.02)

The probabilities of forming a contact are in parentheses.

follows that all top five contacts identified by the two algorithms as well as their rankings are identical apart from a minor reordering of the two contacts Asp23-Lys28 and Ser26Lys28. Moreover, with the exception of Ala21-Val24, the probabilities of contacts are equal within their sampling errors. Finally, visual inspection of the difference contact map ⟨ΔCp(i,j)⟩ = ⟨Cp(i,j)⟩REST − ⟨Cp(i,j)⟩REMD in Figure 8b confirms that the side chain interactions in REST and REMD are in good agreement (all |⟨ΔCp(i,j)⟩| ≲ 0.1). Due to a large number of contacts with small ⟨Cp(i,j)⟩, the relative RMSD comparing REST and REMD contact maps is relatively large (rRMSD = 0.35), but it is reduced to 0.15 if small ⟨Cp(i,j)⟩ < 0.1 max(⟨Cp(i,j)⟩) are excluded. ρ is 1.23, which is close to an ideal value of 1. Compared to Aβ25-35 in water, REST simulations of the Aβ10-40 peptide have an added complexity of sampling peptide−lipid interactions. Moreover, REST simulations keep water and lipids cold at 330 K for all replicas simulated, which may hinder conformational sampling as the peptide is confined to an artificially viscous environment even at high temperatures. Our previous REMD simulations32 have found that Aβ sequence regions S1 and S3 remain bound to the DMPC bilayer surface, whereas hydrophobic regions S4 and partially S2 penetrate into the bilayer hydrophobic core. To investigate the ability of REST to reproduce these previous findings, we plot in Figure 9a the probabilities P(z;i) for amino acid i to occur at the distance z from the bilayer midplane. In addition, superimposed in this plot are the equilibrium positions, ⟨z(i)⟩, of amino acids obtained from REST and REMD simulations. The REST- and REMD-derived positions, ⟨z(s)⟩, averaged over the sequence regions s are presented in Table 4. Together, Figure 9a and Table 4 confirm that the hydrophobic Aβ region S4 in REST simulations penetrates into the bilayer hydrophobic core below the average position of phosphorus atoms zP = 17.35 Å, thus matching the REMD results. As in REMD, region S2 is localized just above the phosphorus atoms. Also in agreement with REMD, hydrophilic S1 and S3 regions in REST simulations remain at the water−lipid interface well above zP. The relative RMSD between the two ⟨z(i)⟩ distributions is just 0.05 indicating minimal differences in ⟨z(i)⟩ between REST and REMD. The corresponding ρ is also small (0.36), further demonstrating that the differences between REST and REMD averages are below the sampling error in REMD. From these data, we deduce that the REST algorithm reproduces well the equilibrium location of Aβ peptides bound to the DMPC bilayer. To check the ability of REST to sample the interactions of amino acids with lipids, we analyzed the contact map ⟨Cl(i,k)⟩,

Figure 9. (a) REST-produced probability distributions P(z;i) for Aβ10-40 amino acids i to occur at the distance z from the bilayer midplane. Black and red thick lines show the average amino acid positions, ⟨z(i)⟩, computed from REST and REMD32 simulations, respectively, which are in excellent agreement. The dashed line marks the average location of the center of mass of phosphorus atoms zP. (b) Number densities of DMPC lipid heavy atoms, nl(r,z), as a function of the distance r to the Aβ10-40 peptide center of mass and the distance z to the bilayer midplane. Vertical black dashed lines separate proximal (r < Rc) and distant (r > Rc) bilayer regions. The bilayer boundary zb(r) is shown by a thick black line. REST- and REMD-produced32 density distributions shown in the right and left halves of the plot are virtually indistinguishable. Both reveal a lipid density depression beneath the center of the Aβ binding footprint at r < Rc.

Table 4. Positions (Å) of Aβ10-40 Sequence Regions along Bilayer Normal algorithm

S1

S2

S3

S4

REST REMD32

24.3 ± 1.2 24.1 ± 1.1

19.0 ± 3.1 17.8 ± 2.8

21.8 ± 1.9 20.7 ± 2.1

16.5 ± 3.4 15.8 ± 3.4

which probes the interactions between amino acid i and the lipid structural group k = G1−G5 (Figure S3). The five most stable contacts derived from REST and REMD are shown in Table 5. Although the ordering of these contacts slightly differs, four out of five top contacts between REST and REMD are identical, and they occur with the probabilities equal within the sampling errors. When computed for the peptide−lipid contact maps ⟨Cl(i,k)⟩, the relative RMSD is found to be 0.20, whereas ρ = 0.7, arguing that the differences between REST and REMD are again small. Hence, REST reproduces Aβ−lipid interactions in agreement with REMD. J


Article


obtained for both distant and proximal lipids. Using REST, we found that the average −⟨SCD⟩ values computed for these two classes of lipids are 0.160 ± 0.001 and 0.125 ± 0.005, respectively. These results almost coincide with the REMDderived values (0.160 ± 0.001 and 0.128 ± 0.002).33 Although Figure 10 shows somewhat diverging −⟨SCD(i)⟩ for carbons i = 4 and 5 in the proximal region, the relative RMSD computed between REST and REMD for the distant and proximal regions are 0.01 and 0.04, emphasizing that, on average, the differences between REST and REMD are negligible. Furthermore, ρ for these two regions are 0.93 and 1.19, respectively. Given the approach of ρ to the value of 1, we conclude that the differences between REST and REMD are about equal to the REMD sampling error. Taken together, we surmise that the use of cold lipids does not appreciably affect Aβ binding to the zwitterionic DMPC bilayer nor the lipid structure. More generally, our data strongly argue that REST and REMD produce nearly identical conformational ensembles of the Aβ−DMPC system at 330 K.

Table 5. Five Most Frequent Contacts between Aβ10-40 and Lipidsa

a

rank

REST

REMD32

1 2 3 4 5

Asn27-G1 (0.75 ± 0.08) Asn27-G2 (0.65 ± 0.15) Lys16-G2 (0.61 ± 0.23) Ala21-G3 (0.60 ± 0.21) Gly38-G1 (0.53 ± 0.07)

Asn27-G1 (0.68 ± 0.07) Lys16-G2 (0.65 ± 0.14) Asn27-G2 (0.58 ± 0.13) Ala21-G3 (0.56 ± 0.18) Ala30-G2 (0.54 ± 0.20)

The probabilities of forming a contact are in parentheses.

Even though our results suggest that REST adequately samples peptide−lipid interactions, the lipid conformations themselves may still be sampled insufficiently. To assess the equilibrium structure of the DMPC bilayer, Figure 9b compares the number densities nl(r,z) of lipid heavy atoms as a function of the distance r to the peptide center of mass and the distance z to the bilayer midplane, which were produced by REST and REMD. This figure shows that REST- and REMD-derived distributions nl(r,z) are nearly identical. More specifically, in agreement with our previous REMD-based findings,32 nl(r,z) computed using REST reveals a depression in lipid density beneath the center of Aβ binding footprint. The extent of the lipid density depression can be quantitatively assessed by defining the bilayer boundary zb(r) (see the Supporting Information). If the boundary zb(r) is computed for the upper and lower leaflets, it affords a straightforward evaluation of the bilayer thickness, D(r), as a distance between upper and lower boundaries. The thinning of the bilayer, ΔD, caused by Aβ binding is then a change in the bilayer thickness between distant and proximal regions; i.e., ΔD = D(r > Rc) − D(r < 6 Å). The resulting ΔD computed for REST and REMD simulations are in excellent (7%) agreement (11.8 ± 3.1 vs 12.7 ± 2.7 Å). As a final probe of lipid conformations in the DMPC bilayer, we computed the lipid carbon−deuterium order parameters −⟨SCD(i)⟩ for each carbon i in the sn-2 fatty acid chains of DMPC lipids (Figure 10). The distribution −⟨SCD(i)⟩ was

4. DISCUSSION In this study we compared the quality and efficiency of REST conformational sampling against the REMD controls. We showed that, compared to REMD, the number of replicas in REST simulations can be reduced four to five times without changes in the temperature range (from 16 to 4 for Aβ25-35 in water or from 40 to 8 for Aβ10-40 binding to the DMPC bilayer). Despite such a dramatic decrease in the number of replicas, the rate of replica exchanges between temperatures was not compromised remaining at about 20−30% due to significant overlaps in enthalpy distributions (Figure 2). Furthermore, m(T) probing distributions of replicas over temperatures approach their optimal theoretical values for Aβ25-35 and Aβ10-40 systems as seen in Figure 3, thus demonstrating nearly ideal replica−temperature mixing in REST. The last conclusion is directly confirmed by a random walk, which replicas perform throughout REST trajectories (Figure S2). Nearly ideal mixing of replicas over temperatures is a necessary condition for exhaustive conformational sampling. To probe the extent of conformational sampling, we considered Ns collected with simulation progress. For both (small and large) simulation systems, Ns collected from all REST replicas in Figure 4 levels off indicating an exhaustion of new states. Although this behavior qualitatively mirrors that of REMD, the final numbers of unique states collected by REST are up to 10-fold smaller than the REMD corresponding numbers. (Importantly, for computing Ns, REST enthalpies scaled according to eq 1 were “unscaled” back to recover wildtype values.) An obvious reason contributing to the dramatic gap in the numbers of unique states produced by REST and REMD is a much smaller number of REST replicas. However, even for the same τsim, REST-producedNs are still much smaller than their REMD counterparts. A reason for this difference is that in REST replicas the bulk of the system (water or water with lipids) remains cold at all temperatures, thus limiting the range of enthalpy values. In this context, it is important that, despite decreased diversity of conformational states in REST, when the collection of new states is restricted to a target temperature (320 K for Aβ25-35 or 330 K for Aβ10-40), the distinctions between REMD and REST disappear or are drastically reduced (insets to Figure 4). This observation indicates that the REST replicas featuring scaled (non-wildtype) energy functions do not significantly degrade the

Figure 10. Comparisons of lipid carbon−deuterium order parameters, −⟨SCD(i)⟩, computed for carbon atoms i in the fatty acid tails sn-2. Data in black and red are obtained from REST and REMD32 simulations, respectively. Continuous and dashed lines represent −⟨SCD(i)⟩ computed for the proximal and distant regions. Vertical bars indicate sampling errors. The figure shows that proximal and distant lipids are sampled in REST and REMD simulations equally despite being kept “cold” in the former. K


Article


(REST) or 94% (REMD) contribution to the wild-type (330 K) partition function. Thus, nearly identical sampling of states at 330 K alone produced by REST and REMD is sufficient to generate conformational ensembles that are in perfect agreement. It should be pointed out that potential utility of REST for simulating lipid bilayers was predicted by Berne and coworkers.22 Our study has confirmed that theoretical expectation. For the Aβ25-35 peptide in water, the REST sampling errors (especially for the turn fraction per residue) are larger than their REMD counterparts implicating much slower convergence of REST simulations than for Aβ10-40. We found that, depending on the specific structural quantity, the number of REST trajectories should be increased up to four times compared to REMD to reach equal sampling errors. Using relative RMSD values, we established that the turn propensities between REST and REMD simulations differ, on average, by 20%, whereas the probabilities of occurrence of backbone dihedral angles deviate by about 37%. The difference in the probabilities of side chain contacts is around 34%, and about 21% difference is observed in the probability distributions of radius of gyration. Although the REST- and REMD-produced Aβ25-35 conformational ensembles are qualitatively similar, the differences between them are clearly larger than those seen for the Aβ10-40 system. Moreover, ρ for all Aβ25-35 quantities examined were larger than 1.7 reaching as much as 2.7 for turn propensities. These data indicate that REMD and REST conformational ensembles are statistically different. Then, a question arises if REST and REMD sampling will ever converge. To answer this question, we consider the relative RMSD, rRMSD, comparing REST turn propensities against full REMD sampling as a function of the number of REST trajectories, M. Because REMD sampling uses all available trajectories (MREMD = 5), rRMSD(M) probes the approach of REST simulations to the REMD sampling. Figure 11 shows

sampling of the wild-type replica, at least in terms of collecting new states. Although this circumstance alone does not imply that REST and REMD conformational ensembles would be identical, it represents a necessary requirement for their consistency. We directly compared the sampling errors in REST and REMD simulations. As expected in both systems, Aβ25-35 and Aβ10-40, these errors decrease with the number of independent trajectories for all structural quantities considered. It is interesting to note that for Aβ peptides binding to the DMPC bilayer the sampling errors in the helix propensities of amino acids or in their positions along the bilayer normal (δH and δz in Figure 5c,d) are similar between REST and REMD. To ascertain that REST and REMD produce consistent Boltzmann weights of conformational states, we analyzed the thermodynamic averages of multiple structural probes for the Aβ10-40 system. Using relative RMSD values, we have determined that the REST- and REMD-produced helical propensities for Aβ10-40 peptides binding to the DMPC bilayer agree, on average, within 10% (rRMSD = 0.10), whereas the differences in the probabilities of forming side chain contacts are about 15% (excluding contacts with negligible probabilities of occurrence). An excellent agreement between REST and REMD is observed in the locations of amino acids along the bilayer normal (within 5%), although somewhat larger differences between REST and REMD (20%) characterize amino acid−lipid interactions. The deviations in the bilayer thinning and carbon−deuterium order parameters of fatty acid tails are 7% and up to 4%, respectively. These results argue that REST reproduces REMD thermodynamic averages extremely well. Several comments about the REST simulations of Aβ10-40 binding to the DMPC bilayer are noteworthy. First, ρ computed for Aβ10-40 and Aβ-lipid interactions is typically smaller than 1 suggesting that the differences between the REST and REMD averages are actually smaller than the REMD sampling error. This outcome can be rationalized as follows. Because REST and REMD trajectories were started with the same initial structures, small ρ values (

Does Replica Exchange with Solute Tempering ... - ACS Publications

Recommend Documents