Balancing Accuracy and Cost of Confinement Simulations by

Apr 27, 2016 - Improvements to the confinement method for the calculation of conformational free energy differences are presented. By taking advantage...
1 downloads 7 Views 2MB Size
Download PDF
Article pubs.acs.org/JCTC

Balancing Accuracy and Cost of Confinement Simulations by Interpolation and Extrapolation of Confinement Energies François Villemot,† Riccardo Capelli,‡,§ Giorgio Colombo,‡ and Arjan van der Vaart*,† †

Department of Chemistry, University of South Florida, 4202 East Fowler Avenue, CHE 205, Tampa, Florida 33620, United States Istituto di Chimica del Riconoscimento Molecolare, Consiglio Nazionale delle Ricerche, Via Mario Bianco 9, 20131 Milano, Italy § Dipartimento di Fisica, Università degli Studi di Milano and INFN, Via Celoria 16, 20133 Milano, Italy ‡

ABSTRACT: Improvements to the confinement method for the calculation of conformational free energy differences are presented. By taking advantage of phase space overlap between simulations at different frequencies, significant gains in accuracy and speed are reached. The optimal frequency spacing for the simulations is obtained from extrapolations of the confinement energy, and relaxation time analysis is used to determine time steps, simulation lengths, and friction coefficients. At postprocessing, interpolation of confinement energies is used to significantly reduce discretization errors in the calculation of conformational free energies. The efficiency of this protocol is illustrated by applications to alanine n-peptides and lactoferricin. For the alanine-n-peptide, errors were reduced between 2- and 10-fold and sampling times between 8- and 67-fold, while for lactoferricin the long sampling times at low frequencies were reduced 10−100-fold.

■

INTRODUCTION Conformational free energy differences are of considerable importance to protein science because they can be used to rationalize and predict the relative stabilities of different protein folds, to elucidate the mechanism of allostery and conformational transitions, and to explain the effect of mutations on these processes. Obtaining conformational free energies by computational methods is no easy task, however.1−8 Most methods require a pathway in configuration space that connects the one state of interest to the other. Pathways with physiological relevance, that is, with high statistical weights and low free energy barriers, are hard to determine due to the enormous dimensionality of configurational space. Moreover, while many methods use a lower-dimensional order parameters to describe pathways, it is generally difficult to predict and verify if these order parameters are sufficient and correct.9,10 Other pathways, like interpolated pathways, suffer from the occurrence of large free energy barriers, which impede sampling.11,12 The difficulties of calculating pathways and associated free energy profiles in configurational space can be circumvented by the use of nonphysical pathways to connect the states of interest. This principle is exemplified by the confinement method,13−19 which connects the states of interest by artificial, uncoupled harmonic oscillator states. In the confinement method, the states of interest are slowly transformed into independent harmonic oscillators using harmonic restraints of increasing strength. Because the free energy cost of applying these restraints can be readily calculated and the free energy of the harmonic oscillator is known, conformational free energy differences can be obtained in a relatively straightforward © XXXX American Chemical Society

manner, especially when using a modified harmonic oscillator state that is invariant to rigid body motions.13,14 This calculation is normally carried out in an implicit solvent, but the method has recently been extended to include explicit solvation.15 By foregoing a path in configurational space, the confinement method provides an efficient and robust way to calculate conformational free energy differences, even for states that are highly dissimilar in structure, but the method has other computational benefits that can be exploited. The free energy of transforming the system to a set of independent harmonic oscillators is obtained through a series of restrained simulations, each with a different strength of the harmonic restraint. While the strengths differ, the center of the restraint is the same in all of these simulations and corresponds to the equilibrium structure of interest. This means that all configurational space accessible to the high restraint simulation is contained in the accessible configurational space of the low constraint simulation. Naturally, the numerical value of this overlap will be very small when there is a large difference in restraints (because the amount of accessible configurational space shrinks as the restraint increases in strength); however, the overlap can be particularly large for “neighboring” simulations for which the restraint force constants are closest in value. The overlap in configurational space is closely associated with the overlap of energy distributions, which is crucial for the accurate estimate of thermodynamic properties.20 Received: December 14, 2015

A

DOI: 10.1021/acs.jctc.5b01183 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

between successive points, would therefore be a bad choice. Because this relation becomes linear in logarithmic space, the integration is carried out using linear interpolation in logarithmic space instead;21 in the following, we will refer to the integrand as

Here we exploit this overlap in order to maximize the efficiency and minimize the statistical error of confinement simulations. We will show that one can use the overlap in configurational space to accurately predict confinement energies at unsampled strengths and that this interpolation significantly decreases the error in calculated free energies. We will also show that instead of sampling at given intervals, one can use the overlap to predict at which restraining strength to sample next for simultaneously optimal errors and costs. Finally, by coupling these interpolations and extrapolations to relaxation time analyses, we will introduce a robust protocol for optimal confinement simulations, which is illustrated by applications to the alanine n-peptide and lactoferricin.

Iν ≡

■

∫0

1

NDOF log(βhν) − β

∫0

Uconf ν2

dζ

(4)

The averages X0)⟩Hζ(ν2) are obtained from independent simulations that generate representative configurations at the given restraint frequencies. Conformational Free Energies. By removing rigid body motions in the confinement procedure,14 GΩ of eqs 1 and 2 represents a vibrational configurational free energy that lacks free energy contributions from the overall translational and rotational motions. The translational free energy can be obtained from the partition function of the ideal gas, while the rotational component can be obtained from the partition function of the rigid rotor.22,23 In calculating conformational free energy differences (ΔG) between two states, the translational components will cancel, but rotational components generally will not, because of differences in the moments of inertia. Here, all reported ΔG values include the rotational component, while ΔGΩ is used to denote a difference in vibrational free energy. Reweighting and Interpolation. Because the system is confined to the vicinity of the same reference structure X0 in each simulation, there is large spatial overlap between these sets of configurations. This means that the configurations obtained at a given frequency can be used to estimate ensemble averages at a different frequency. Consider Ni configurations obtained from a simulation with restraint frequency νi. The ensemble average of observable A at frequency νj is given by

(1)

ν2

NDOF (2π )2 βM

⟨ρ2m(X,

Here, β is the inverse temperature (β = 1/kbT), h is Planck’s constant, M is the total mass of the system, ρm (X, X0) is the mass-weighted root-mean-square distance (rmsd) of sampled configuration X with reference structure X0, E(X0) is the potential energy of the reference structure, and ⟨.⟩ denotes an ensemble average. The Hamiltonian Hλ(λ) depends on the parameter λ, which is used to switch from an unrestrained system (λ = 0) to a system with harmonic restraints of frequency ν (λ = 1). This frequency must be chosen high enough so that at λ = 1 virtually all of the system’s energy is due to the restraints. At that point, the system can be considered to be a set of noninteracting harmonic oscillators. For convenience, the integration in eq 1 can be transformed by the change of variable ζ = λν2 and Hζ(X;ζ) = Hλ(X;λν2), so that the integration is done in frequency space. In addition, the mass-weighted rmsd can be expressed as a function of the confinement energy Uconf = 2Mπ2ν2⟨ρ2m(X, X0)⟩Hζ(ν2): GΩ = E(X 0) +

(3)

ν 2⟨ρm2 (X , X 0)⟩Hζ(ν 2) =

NDOF log(βhν) − 2(πν)2 M β

⟨ρm2 (X , X 0)⟩Hλ(λ) dλ

ν2

In practice, simulations are carried out with increasing values of ν until the kinetic energy of the system (NDOFkbT/2) equals the confinement energy Uconf, as expected for a purely harmonic system, or equivalently:

METHODS Confinement Method. The confinement method13−19 aims to compute the free energy of macrostate Ω by transforming the system of interest with 3N degrees of freedom into a set of 3N noninteracting harmonic oscillators (HOs). This transformation is performed in a series of simulations by adding harmonic restraints of increasing strength to the physical potential. The restraints are centered at the positions of a reference structure X0 belonging to Ω, and the simulations are performed until the system can be considered purely harmonic. A recent improvement,13,14 also used here, removes the rigid-body motions by performing a mass-weighted best-fit alignment of the system onto X0 at each simulation step. Removing the sampling of the translational and rotational degrees of freedom significantly speeds up convergence, especially at low restraint strengths. The total number of degrees of freedom (NDOF) then becomes 3N − 5 for systems with linear geometries and 3N − 6 otherwise. The vibrational free energy of the Ω macrostate can subsequently be computed by thermodynamic integration:14 GΩ = E(X 0) +

Uconf

⟨A⟩νj = ⟨A e−β(Uj − Ui)⟩νi e−β ΔFij

(5)

Here ΔFij = Fi − Fj is the free energy difference between the biased states at νi and νj; we use the symbol F to distinguish it from the configurational free energy of the unbiased state (G) of eqs 1 and 2. Uj and Ui are the potential energy values corresponding to frequencies νj and νi, respectively. Because only the restraint differs between these potentials, ⎡⎛ ν ⎞ 2 ⎤ j Uj(Xk) − Ui(Xk) = 2Mπ 2ρm2 (Xk , X 0)νi2⎢⎜ ⎟ − 1⎥ ⎢⎣⎝ νi ⎠ ⎥⎦

(6)

the accuracy of the ensemble average obtained by reweighting (eq 5 and 6) diminishes with increasing |νj − νi|. A more accurate estimation can be obtained by using additional statistics. This can be done by combining configurations obtained from all simulations at all frequencies. To combine samples from these multiple simulations, an estimation of the free energy difference between the states is needed, which can be obtained from the multistate Bennett acceptance ratio estimator (MBAR).24

(2)

For a harmonic oscillator, ⟨ρ2(X, X0)⟩ ∝ ν−2. A trapezoidal rule for numerical integration of eq 2, which interpolates linearly B

DOI: 10.1021/acs.jctc.5b01183 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

conformations of the alanine dipeptide are C7ax and C7eq, which for the force-field and implicit solvent method used, correspond to backbone dihedral angles of (ϕ, ψ) = (61.4, −71.4) and (−78.0, 138.7) degrees, respectively. The C7ax and C7eq conformations were used as the reference structures for the alanine dipeptide. Larger alanine n-peptide systems behave nearly like independently linked alanine dipeptides when the peptides are in the C7ax or C7eq states.12,29 For instance, the energy minima (C7ax, C7ax) and (C7eq, C7eq) of the alanine tripeptide correspond to (ϕ1, ψ1, ϕ2, ψ2) = (61.1, −72.1, 59.6, −71.6) and (−77.4, 137.7, −76.7, 137.8) degrees, respectively. For all n ≥ 3 alanine systems, the confinement reference structures were obtained by setting all the (ϕ,ψ) dihedral angles to the values of C7ax and C7eq of the alanine dipeptide and performing an energy minimization. In the following, these configurations are simply named C7ax and C7eq, independently of the number of dihedral angles. The moments of inertia were obtained for the reference structures. The corresponding contribution to the free energy difference between C7ax and C7eq equals

With this reweighting, we can interpolate averages between simulated frequencies, thereby increasing the accuracy of the thermodynamic integration of eq 2 in a cost effective way. We also use it to extrapolate the value of the confinement energy for frequencies higher than the highest simulated thus far. As described below, this extrapolation is used to assess the frequency of the next simulation such that the overall cost and accuracy of the procedure is optimized. Finally, the reweighting also increases the accuracy at the simulated frequencies by mixing in configurations from the other simulations in the computation of the confinement energy. Extrapolation. To properly estimate the error on the extrapolated value, the following setup was used: we start from a set of simulations performed with different restraints, up to a frequency νmax. We subsequently employ the extrapolation of eqs 4 and 5 to estimate the confinement energy for a set of unsampled frequencies νi > νmax. The computational cost of this extrapolation is low (much lower than the actual sampling) and nearly independent of the number of unsampled frequencies. We chose this direction because, in extrapolating toward higher frequencies, phase space is compressed. This means that all relevant areas of space for the higher frequency restraint were sampled in the lower frequency simulation (but insufficiently). If we were to choose the other direction, that is, sampling at the higher frequency followed by extrapolation to the lower frequency, certain regions of space important for the low frequency restraint would be left unsampled. After the extrapolation, a confinement simulation is performed for each of the new frequencies in order to obtain the actual value of the confinement energy, and these calculated values are compared with those obtained from the extrapolation. The ratio between the extrapolated and actual confinement energy is a measure of the error of the extrapolation. We express this ratio as a function of the free energy difference ΔF between the simulations at νi and νmax. ΔF can be obtained from eqs 4 and 5, or, if simulations at multiple frequencies are used, from MBAR. The latter approach would yield somewhat more accurate extrapolations because more data is used. However, here we used data from only one simulation and the former approach in order to base all comparisons on the same amount of data. For a given νmax, ΔF increases with |νmax − νi| and represents a meaningful quantity that can be compared across systems of different sizes. Correlation Times. The efficiency and accuracy can be improved further by considering the correlation time of the system. This correlation time is affected by the addition of harmonic restraints, especially at high frequencies, when the confinement energy accounts for a large portion of the total potential energy. In addition, these restraints limit the configurations accessible to the system to the ones close to the reference structure X0, and the phase space to sample gets smaller as the frequency gets higher. To attain comparable sampling for each frequency, different sampling times are therefore needed, which can be estimated from the correlation time. These were estimated by block-averaging the confinement energies25 and also by calculating the autocorrelation function of the confinement energy. Alanine n-Peptide Setup. We performed confinement simulations of capped alanine n-peptides (n = 2, 4, 6, 8, and 10), with the general formula CH3CO-Alan−1-NHCH3. These simulations were carried out with the CHARMM program,26 using the CHARMM polar hydrogen parameter set param1927 and the ACE implicit solvent model.28 The two lowest-energy

ΔGrotation =

⎛ ⎞ k bT ⎜ IC7ax ⎟ ln⎜ ⎟ 2 ⎝ IC7eq ⎠

where IC7ax and IC7eq represent the products of the three principal moments of inertia of C7ax and C7eq, respectively.22,23 All confinement simulations were performed using Langevin dynamics at 300 K, with friction coefficients of 1, 5, 10, or 20 ps−1 (see Results). The time step had a maximum value of 1 fs and was adjusted depending upon the restraint frequency. It was chosen so there are at least 80 time steps per harmonic oscillator period, resulting in smaller time steps for higher frequencies. Different time steps were tested, which showed that at least 40 steps per period are required to obtain accurate estimation of the confinement energy. A conservative value of 80 steps/period was then chosen. SHAKE30 was not used in the simulations. To further restrict sampling to the state of interest, especially at the lowest frequencies, we also added flat-bottom dihedral restraints. These were centered on the energyminimized values, with a force constant of 10 kcal/mol/rad2 and a width of 2.5°. This value was chosen so that the states of interest are the same as the ones defined in the umbrella sampling. Interpolation of the confinement energy was done for 10 frequencies, equally spaced in log-space, between consecutive simulations. Adding more points did not change the final free energy difference or the error bars. For comparison, free energy differences between the C7ax and C7eq, conformations were also obtained from one-dimensional umbrella sampling simulations.31 For the alanine n-peptide, the transformation from C7eq to C7ax involves (n − 1) ϕ and (n − 1) ψ dihedral transformations. These angles were treated as reaction coordinates and changed one at a time (in the order ϕ1,ψ1,⇌,ψn−1) while keeping the others constant. For each dihedral angle, 50 equally spaced umbrella windows were used, with a force constant of 150 kcal/mol/rad2 and a simulation time of 100 ns per window. To maintain the trans peptide configuration, flat-bottom dihedral restraining potentials were used for the ω backbone dihedral angles with a force constant of 10 kcal/mol/rad2 and a width of 90°. Simulations were performed with Langevin dynamics at 300 K, using a 1 fs time step, no SHAKE,30 param19,27 and ACE.28 Potentials of mean C

DOI: 10.1021/acs.jctc.5b01183 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation force (PMF) were obtained from MBAR;24 all free energies are reported in kcal/mol. Lactoferricin Setup. Bovine lactoferricin is 25-residue peptide cleaved from lactoferrin with antimicrobial properties.32 In lactoferrin, the sequence is folded into an α-helix followed by a β-strand, while the cleaved peptide adopts a β-hairpin fold; the peptide contains one disulfide bond (Figure 1).33,34 No

solvent model,37 Langevin dynamics, and no SHAKE.30 A friction coefficient of 1 ps−1 was used for simulations with a frequency lower than 2 ps−1 and a friction coefficient of 20 ps−1 for frequencies above. Interpolation of the confinement energy was done for 10 equally log-spaced frequencies between consecutive simulations. After an energy minimization, each conformation of the peptide was heated and equilibrated at 300 K. The reference structures used in the confinement simulations were obtained from rmsd-based clustering with a cut off of 3.5 Å of a 25 ns unrestrained trajectory. All simulations were conducted with the CHARMM program.26

■

RESULTS Alanine n-Peptide. Figure 2 shows the ratio between the extrapolated and actual confinement energies as a function of ΔF for multiple values of νmax. Curves for all alanine systems are provided; a value of one indicates that the extrapolation perfectly predicted the confinement energy. As expected, deviations from one strongly increased with the free energy difference and the ratio was very close to one for small ΔF. A notable feature is that the extrapolation stayed accurate for larger free energy differences as νmax increases. In other words, as the system becomes more harmonic, it becomes easier to predict the result of a new simulation. This is due to the fact that at larger frequencies, the harmonic restraints represent a larger portion of the total energy and the configurational space is compactly distributed around X0 in a predictable manner. In addition, we observed that the extrapolation is more accurate as the size of the system increases. Larger systems have narrower energy distributions, so that the weights in eq 5 are closer to one another. More configurations will therefore contribute significantly to the ensemble average at another frequency, thus lowering the error. This is a particularly encouraging feature which will facilitate the application of the confinement method to larger and more complex systems. The information on Figure 2 can be used to extract the maximum value of ΔF for which the extrapolation error is below a desired threshold. We chose 5%, a fairly conservative

Figure 1. Structure of lactoferricin. Solution structure on right, structure of the lactoferricin sequence within the lactoferrin protein on left, and disulfide bonds shown as ball and stick.

spontaneous conformational transitions were observed in long unbiased MD simulations.35 Because of its size and the complexity of the transition, lactoferricin is a good test system for the confinement method and representative of the more challenging biological systems that are the ultimate target for the method. The α+β conformation was obtained from residues 17−41 of lactoferrin (PDB: 1BLF34), while the β-hairpin conformation was taken from lactoferricin in solution (PDB: 1LFC33). We used the CHARMM36 force field36 with the GBMV implicit

Figure 2. Ratio between the confinement energy computed from extrapolation and from an actual simulation for alanine 2-peptide (A), 3-peptide (B), 4-peptide (C), 6-peptide (D), 8-peptide (E), and 10-peptide (F) . Simulations with different restraints up to a frequency νmax were used to extrapolate the confinement energy at higher frequencies. These extrapolated states have a higher free energy than the state corresponding to νmax, with a difference of ΔF. D

DOI: 10.1021/acs.jctc.5b01183 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

the extrapolation because more data is available for the interpolation (one extra set of simulations). Furthermore, the error of interpolation can be reduced further by taking into account all simulated data at all simulated frequencies. As shown below, this interpolation significantly reduced the overall error in calculating the configurational free energies. Thus, we propose to exploit ΔFext as guideline for selecting the frequency spacing of the simulations. The goal of this procedure is to pick the maximum spacing at which high quality interpolations remain feasible, thereby obtaining high accuracy at minimal computational costs. If we have a set of simulations up to frequency νmax, the next frequency of simulation will be picked such that its free energy difference with the νmax simulation is ΔFext. Figure 4 shows that interpolation can be performed to obtain confinement energies at nonsimulated frequencies. In Figure 4A, the value of Iν (eq 3) is shown in black for the alanine dipeptide at simulated frequencies of 0.021 and 6.6 ps−1. The free energy difference between these simulations was 4.2 kcal/ mol. The black line represents the value of Iν that would vary linearly with the logarithm of the frequency, which is the assumption made when performing the integration of eq 2 in log space18 and also the analytical solution for harmonic oscillators. The red curve corresponds to interpolated values using MBAR. Additional simulations at intermediate frequencies confirm the accuracy of the interpolation. The simulated values (blue dots) show that Iν does not follow a straight line in log space but falls on the interpolated curve instead. The observed nonlog−linear behavior is expected for this frequency range because the system is far from being purely harmonic. In fact, the harmonic terms contribute only 24% of the total energy at a frequency of 6.6 ps−1. The interpolation accurately reproduced the observed behavior, which demonstrates that meaningful information about the system can be obtained through interpolation. It also shows how the interpolation can greatly increase the efficiency of the method: while all simulations (represented by black and blue symbols) would be needed to accurately compute the free energy difference over that frequency range, just two initial simulations (black point) are sufficient if the interpolation is used. Greater

value for this error, which corresponds to an extrapolated/ actual confinement energy ratio of 0.95 or 1.05. In the following, we will refer to this spacing as ΔFext. Figure 3 shows

Figure 3. Free energy difference corresponding to an error of 5% on the confinement energy obtained by extrapolation. This free energy difference is between the simulated system having the highest harmonic restraints and the one corresponding to the extrapolated frequency.

ΔFext as a function of frequency for the alanine n-peptides. While the curves are bumpy (due to the fact that the ratios switched between 0.95 and 1.05, discretization of ν, and finite sampling), the graph shows several clear trends: consistent with the results of Figure 2, ΔFext increased with both frequency and system size. The physical relevance of ΔFext is the following. When the free energy difference between the sampled system at νmax and the unsampled system at higher frequency is ΔFext, there is sufficient overlap in distribution functions to estimate, within some preselected error bound (here 5%), the confinement energy at the unsampled frequency from simulated data at νmax. This means that after sampling at both frequencies, there will be sufficient overlap in distribution functions to accurately calculate confinement energies at frequencies in between. The accuracy of this interpolation will be higher than the accuracy of

Figure 4. Interpolation of the confinement energy for alanine 2-peptide for different restraint frequency ranges: (A) low frequencies (0.02−6.6 ps−1), (B) intermediate frequencies (6.6−18.4 ps−1), and (C) high frequencies (26.6−38.8 ps−1), as well as (D) for alanine 10-peptide at low frequencies (0.02−0.37 ps−1). The value of the confinement energy is interpolated (red) for frequencies between two initial simulations (black). The thickness of the red line represents the error bar. Additional simulations (blue) are added thereafter to estimate the accuracy of the interpolation. E

DOI: 10.1021/acs.jctc.5b01183 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 5. Correlation times for (A) alanine 2-peptide and (B) alanine 10-peptide for different restraint frequencies, calculated from the autocorrelation function of the confinement energy (triangle symbols) and from block analysis (circle symbols). Multiple friction coefficients (from 1 to 20 ps−1) were used as parameters for the Langevin thermostat.

Table 1. Free Energy Differences between the C7eq and C7ax Conformations of the Alanine n-Peptide in kcal/mola strategy 1b system Ala2 Ala3 Ala4 Ala6 Ala8 Ala10

ΔGUS −3.42 −7.01 −9.38 −16.25 −22.6 −28.98

± ± ± ± ± ±

ΔGHom 0.06 0.09 0.21 0.25 0.31 0.35

−3.57 −6.84 −10.03 −17.21 −23.41 −29.87

± ± ± ± ± ±

0.27 0.35 0.39 0.48 0.46 0.33

costf (%)

ΔG −3.25 −6.91 −9.72 −16.96 −23.47 −30.94

± ± ± ± ± ±

2c

0.06 0.05 0.05 0.05 0.05 0.05

2.57 3.54 4.94 6.72 9.24 12.28

costf (%)

ΔG −3.15 −7.03 −10.15 −16.75 −23.68 −30.69

± ± ± ± ± ±

3d

0.06 0.06 0.08 0.09 0.09 0.09

2.41 2.20 2.33 2.88 3.27 3.67

costf (%)

ΔG −3.77 −6.99 −10.38 −17.20 −23.85 −29.60

± ± ± ± ± ±

4e

0.07 0.07 0.07 0.06 0.06 0.16

2.19 2.29 3.12 4.46 5.50 6.87

costf (%)

ΔG −3.60 −6.67 −10.10 −17.47 −23.74 −30.30

± ± ± ± ± ±

0.11 0.07 0.08 0.10 0.09 0.09

1.50 1.98 2.86 3.67 4.47 5.39

a ΔGUS corresponds to the free energy difference calculated from one-dimensional umbrella sampling. ΔGHom corresponds the free energy obtained from confinement at frequencies that are equally spaced in log space. A total of 17 of these were performed per conformation, with a constant simulation time of 20 ns per simulation, according to the setup of ref 14. In the strategies 1−4, the simulations are spaced in free energy space according to a given relation (see text). All free energies are in kcal/mol. bΔF = 5 kcal/mol. cΔF from system-dependent best fits of to ΔFext (Figure 3). These fits are: ΔFAla2 = 3.19 + 0.54 log(ν/ν1), ΔFAla3 = 1.44 + 1.02 log(ν/ν1), ΔFAla4 = 0.95 + 1.27 log(ν/ν1); ΔFAla6 = 3.36 + 1.12 log(ν/ν1), ΔFAla8 = 3.27 + 1.38 log(ν/ν1), ΔFAla10 = 4.87 + 1.55 log(ν/ν1). dΔFext = 3 + 0.388 ln(ν/ν1), shown by the lower dotted line in Figure 3. eΔFext = 5 + 0.485 ln(ν/ν1), shown by the upper dotted line in Figure 3. fTotal computational cost as a percentage of the total computational cost when using homogeneous spacing in log frequency.

while for the decapaptide, higher friction coefficients led to higher correlation times in this frequency range. This is likely due to the more complex landscape of the decapeptide, which has subbasins; visiting the various subbasins is hindered by large friction terms. Near a frequency of 0.2 ps−1, the correlation times dropped significantly for all systems. At this frequency, Uconf represents between 2 and 4% of the kinetic energy. Apparently, this energy is sufficient to limit the system to one subbasin, which explains the precipitous decline in correlation time. At high frequencies, low correlation times were observed, inversely proportional to the friction coefficient. We checked that the average values of the confinement energies were not affected by the friction coefficient so that higher friction coefficients indeed led to faster sampling. The simulation time needed to obtain a given number of independent measurements is proportional to the correlation time. Figure 5 shows that in order to obtain uniform sampling across the frequencies, much smaller simulation times are needed at higher frequencies. In addition, by increasing the friction coefficient at high frequency, the simulation time can be reduced, thereby further increasing the efficiency of the calculation.

accuracy can also be achieved by reducing the discretization error arising from the frequency spacing, but this comes at additional computational costs. Figure 4B,C illustrate how the interpolation for the alanine dipeptide performs at higher frequencies. The free energy differences between these two frequencies are comparable to the one corresponding to Figure 4A (5.9 and 4.2 kcal/mol versus 4.2 kcal/mol). Again, the interpolation correctly estimated Iν for frequencies that were not simulated. At higher frequencies, Iν varied more linearly with the log of the frequency, as expected for more harmonic system. The same behavior was observed for the other alanine systems, such as alanine 10-peptide (Figure 4D). Figure 5 shows the correlation time of the confinement energy for the alanine dipeptide and decapeptide as a function of the restraint frequency. The correlation times were calculated by block-averaging25 (indicated by circles) and from the autocorrelation function of the confinement energy (triangles). The two methods gave similar results, which indicates that the correlation time could be properly estimated. For the dipeptide, the correlation time was similar for all frequencies