(PDF) Implementation of the Forward–Reverse Method for Calculating the Potential of Mean Force Using a Dynamic Restraining Proto

Article pubs.acs.org/JPCB Implementation of the Forward−Reverse Method for Calculating the Potential of Mean Force Using a Dynamic Restraining Protocol Mostafa Nategholeslam,*,† C. G. Gray,*,‡ and Bruno Tomberli*,§ † Department of Physics and Biophysics Interdepartmental Group, University of Guelph, Guelph, Ontario, Canada ‡ Guelph-Waterloo Physics Institute and Department of Physics, University of Guelph, Guelph, Ontario, Canada § Department of Physics, Capilano University, North Vancouver, British Columbia, Canada ABSTRACT: We present a new sampling and analysis scheme for calculating the potential of mean force (PMF) of systems studied by steered molecular dynamics simulations. This scheme, which we call the bin-passing method, is based on the forward−reverse (FR) method (due to I. Kosztin and co-workers, Kosztin et al. J. Chem. Phys. 2006, 124(6), 064106) and arguments based on the second law of thermodynamics. Applying the bin-passing method results in enhanced sampling, better separation of the reversible and irreversible work distributions, and faster convergence to the underlying PMF of the system under study. Post-simulation analysis is performed using a purpose-built software that we have made publicly available at https://github.com/1particle/bin-passing_analyzer under the terms of the GNU General Public License (version 3). Three examples are provided, for systems of varying sizes and complexities, to demonstrate the efficiency of this method and the quality of the results: for the dissociation PMF of NaCl in water, the bin-passing method obtains PMFs in excellent agreement with that obtained for the same system and using the same force-field through static (equilibrium) methods. The bin-passing method gives a very symmetric PMF for passage of a single water molecule through a DPPC bilayer, and the resultant PMF leads to permeability values in better agreement with experiments than those obtained through previous simulation studies. Finally, we consider the interaction of the antimicrobial peptide HHC-36 with two model membranes and employ the bin-passing method to obtain the PMFs for peptide adsorption to the membranes. The characteristics of these PMFs are consistent with the qualities established for the HHC-36 peptide through in vivo and in vitro experiments, as a non-toxic strong antimicrobial agent. ■ INTRODUCTION The advent of Jarzynski’s equality1,2 and, shortly after that, converged PMF from an equilibrium simulation for the same system,33,34 and thus, some of the anticipated merits of using a Crooks’ theorem3 have led to numerous studies on the NEW theorem are lost. In cases where faster steering speeds are applicability of these theorems for extracting free energy used, higher numbers of trajectories are needed to achieve profiles of many-body systems from the external non- convergence, to the extent that it is often computationally more equilibrium work, both in experimental settings4−6 and convenient to stay as close to the equilibrium limit as possible. simulations (e.g., see refs 7−10). These theorems and their Despite being correct in principle (when sampling over implications have also been the subject of considerable infinitely many paths), Jarzynski’s equality was discovered soon theoretical investigation. Several groups have also employed after its advent to suffer from slow convergence, due to its Jarzynski’s equality to obtain free energy profiles in numerous reliance on very small, zero, and even negative values of biomolecular systems.11−31 dissipative work, all of which are known to occur rarely, from Despite all the insight and utility provided by these works, the second law of thermodynamics.35,36 Various authors the original hopes raised by non-equilibrium work (NEW) (including Jarzynski himself) have suggested that using bi- theorems and their extensions, i.e., to provide fast-converging directional work samplings might partially remedy that methods for extracting the potential of mean force (PMF) are difficulty.35,37 The best-known general bidirectional work not entirely fulfilled for large thermodynamic systems, e.g., theorem is that of Crooks3 biomolecular systems. In particular, while a small ensemble of steered molecular dynamics (SMD)7,9,32 trajectories is often ⟨f ( +W )⟩F generated and used to obtain the PMF, thorough analysis of the e−β ΔG = ⟨f ( −W )e−βW ⟩R (1) convergence of the PMF and proper error analysis is often absent from these works. Furthermore, the simulation time invested to achieve the PMF from non-equilibrium SMD Received: May 20, 2014 simulations using a direct application of Jarzynski’s equality, is Revised: November 4, 2014 comparable, if not longer than the time required to achieve a Published: November 5, 2014 © 2014 American Chemical Society 14203 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B ■ Article where ΔG is change in free energy between the two BIN-CROSSING VERSUS BIN-PASSING SCHEMES macrostates under study, and ⟨...⟩F and ⟨...⟩R denote averages Over the course of an SMD simulation, time is discretized into taken over ensembles of paths in the forward and reverse very short time-steps, typically on the order of a femtosecond. directions, respectively. This theorem by itself is not directly We start at xA at t = 0, and proceed xA → xB with a constant applicable to extract PMF profiles, as it contains an arbitrary velocity of v = (xB − xA)/T, where T is the simulation time we function ( f) of work W. One needs to make an explicit choice invest in each realization of the xA → xB process. The forward for that function before being able to implement Crooks’ and the reverse drift speeds need to be the same for eq 2 to be theorem, and some choices may lead to better convergence applicable, so the system will also be steered with the speed v = than others. The forward−reverse (FR) method,33 due to (xA − xB)/T along the reverse xA ← xB path. The specifics of Kosztin and co-workers, is based on Crooks’ fluctuation the choice of the spring constant, total simulation time, and the theorem (the FR method can also be derived from the steering speed, v, are not the subjects of our direct concern Jarzynski equality, within the same assumptions used to derive here. Rather, we wish to revisit eq 2 with the notions of it from Crooks’ fluctuation theorem) and arrives at a reversible and dissipative works in mind, and also with some remarkably simple result: for a system obeying Brownian insight obtained from the second law of thermodynamics. To dynamics, when we steer the system with a harmonic potential simplify the discussion, we formally define the notion of bin- of sufficient stiffness to proceed with uniform speed from initial crossing, which we will later contrast to the notion of bin- to final macrostates (corresponding to the initial and final passing. values of the reaction coordinate), the free energy difference Bin-crossing. During an SMD simulation (e.g., an FR run), between those states is given by half the difference of the one bin-crossing occurs when the value of the chosen reaction averages of forward and reverse works, i.e., coordinate of the system (e.g., the radial distance between two ions) is steered to enter one boundary of a bin, pass through ⟨W ⟩F − ⟨W ⟩R the bin, and eventually emerge from the opposite boundary of ΔG = the bin. This marks the completion of one bin-crossing for that 2 (2) specific bin. A schematic view of this process is depicted in Like other non-equilibrium work theorems, the convergence Figure 1. The elementary scheme for recording work values rate of the RHS (right-hand side) to LHS in eq 2 differs from system to system, based on the size of the system, its internal dynamics, and the thermodynamic parameters of the simulation. Equation 2 has been used to extract PMFs for several test systems of various levels of complexity,33,38−41 and its convergence rate has been shown to exceed that of Jarzynski’s equality, when tested on the same systems. While eq 2 already shows a superior convergence rate when trivially implemented, we argue below that much more can be learned from this simple equation and demonstrate how that knowledge can be utilized toward significantly improving the convergence by devising a more careful analysis scheme. In the elementary data-collection and analysis protocol for using eq 2 in a steered molecular dynamics (SMD) simulation, Figure 1. Progress of the system viewed at time-step resolution, for a it is commonplace to divide the reaction path into the desired system that is steered using a harmonic restraining potential in the A number of bins.7,9,33 Higher number of bins corresponds with a → B direction over the [0, T] time interval. In the enlarged view on higher resolution for the PMF to be calculated. Over the course the right, the horizontal dotted lines represent successive time-steps, each of length δt. The system evolves discretely, one time-step at a of the simulation, we record the work done on the system in time. Although the overall progress of the system is from A to B there each bin as it proceeds from initial state A to the final state B are many time-steps during which the progress of the system is in the (and back), corresponding to initial xA and final xB values of the A ← B direction, mostly as a result of thermal fluctuations. reaction coordinate. In traveling from A to B, the value of the reaction coordinate will traverse each bin at least one time. If we repeat this process a large number of times in the forward during an SMD run is based on bin-crossing: we (re)start (F) and reverse (R) directions, we obtain the average works on accumulating the external work done on the system (in either the RHS of eq 2. The change in free energy across each bin forward or reverse direction) as a bin is entered, accumulate all calculated this way can then be summed to obtain the complete the work in subsequent time-steps spent within the boundaries PMF from A to B. We call this elementary scheme the bin- of that bin, up to the time-step when the opposite boundary is crossing method and compare its accuracy with the bin-passing exited. Although this might seem to be the most intuitive way method, which we present in the next section. We then detail a to measure work for calculating the work averages in eq 2, a method for estimating the uncertainties in the PMFs obtained closer look suggests otherwise. The derivation of eq 2 gives also through the bin-passing method, followed by results for the bin- a relation for the average dissipative work (in either direction) passing method applied to three different systems, in increasing ⟨Wd⟩33 order of complexity. These results allow us to make ⟨W ⟩F + ⟨W ⟩R comparisons between the convergence and accuracy of the ⟨Wd⟩ = 2 (3) bin-passing method compared to the bin-crossing method and comment on each method’s applicability in modeling The meaning of eqs 2 and 3 is quite simple: the average bimolecular processes. forward work done on the system, in general consists of a 14204 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article reversible part (equal to the free energy difference) and a dissipative part ⟨W ⟩F = WF,reversible + ⟨Wd⟩ (4) A similar relation can be written for ⟨W⟩R ⟨W ⟩R = WR,reversible + ⟨Wd⟩ (5) By definition, we have WF,reversible = −WR,reversible = ΔG, and for most systems, ⟨Wd⟩F = ⟨Wd⟩R. Exceptions include asymmetric Brownian objects with a constraint on their orientation along the reaction path. Such exceptions aside, for Figure 2. A sample section of the trajectory traversed in an actual most Brownian systems, one can show33 that the average SMD simulation. The track shows the normal distance (x) between dissipative work ⟨Wd⟩ done on the system in forward and the center of mass of an HHC-36 peptide molecule (sequence reverse directions is the same, when pulls in both direction are KRWWKWWRR)43 from the center of mass of a POPE−POPG done with constant and equal drift speeds, using a harmonic membrane patch. The HHC-36 molecule is steered to go toward the membrane during the shown interval, but it can be seen that during steering restraint of large enough stiffness. The second law of several time-steps the progress of the molecule is toward larger x- thermodynamics asserts further for the average dissipative work values (marked red on the graph, and denoted further by colored bars that ⟨Wd⟩ ≥ 0, with the equality holding only for infinitely slow beneath them). A harmonic potential with k = 8000 kcal/(mol·Å2) is pulls (v → 0), during which system always remains close to used to move the peptide, with another harmonic potential of k = equilibrium. This, together with the discrete nature of SMD 25000 kcal/(mol·Å2) keeping the center of mass of the membrane in pulls, causes most of the time-steps to have a positive value of place. Each segment of the trajectory shown here corresponds to dissipative work, although there is nothing in the second law displacement of the peptide (relative to the membrane) during one against any single time-step having zero or negative Wd, even time-step of duration 2 fs. when the system is far from equilibrium.42 From the enlarged section at the right side of Figure 1, it is clear that when system adding to the overall value of the dissipative work in the is steered to proceed in the forward (reverse) direction with a forward direction. stiff spring, it very often happens to go in the reverse (forward) In signal-noise language, and regarding our sampling of the direction. This can be due to thermal fluctuations in the system reversible work as signal and that of the dissipative work as (e.g., impulses from water molecules in an aqueous medium) or noise, the work recorded during each time-step has a a result of corrective steerings when the system has jumped contribution to signal (reversible work) and some to the ahead (overshot) during an earlier time-step. To understand noise (dissipative work). Counting the work from red segments the latter effect, consider that steering is performed using a as forward work cancels some of the signal from blue segments, spring of constant stiffness, which may result in the system while still contributing to the noise, and neither of these is proceeding further than its target position by the end of some desirable. In this same language, we can say that on the RHS of time-steps. This should be corrected in subsequent time-steps, eq 2, the average noise (dissipative work) from ⟨W⟩F cancels according to the predefined trajectory we wish the system to that in ⟨W⟩R, while the signal (reversible work) in ⟨W⟩F sums follow with the negative of that in ⟨W⟩R, which is in turn the negative x − xA of the signal value in ⟨W⟩F, and then dividing by 2 we obtain xtarget(t ) = xA + v ·i·δt , v = B the average value of the signal. T (6) Adopting a scheme where red-segment time-steps are with a similar equation holding in the reverse direction. Here i excluded from our samplings of forward work, does not mean is the time-step counter and δt is the length of each time-step, that we discard them. Rather, they can be used in our sampling such that t = i·δt. of the reverse work ⟨W⟩R, as they contain the correct value of An example of this fluctuation phenomenon in an actual the A ← B reversible and dissipative works. The rationale SMD run is shown in Figure 2, where the system is restrained discussed above holds also in the reverse direction: only time- (or targeted) to proceed to smaller x values (blue segments), steps during which the system actually proceeds in the reverse but every now and then it moves in the opposite direction (red direction should be used in our estimation of ⟨W⟩R, and these segments, denoted also by colored bars beneath them). can occur both during the (overall) forward and reverse portion Our major assertion here is that the external work done on of the SMD run. the system over the red segments in Figure 2 should not be Although these assertions are intuitively appealing, imple- counted as forward work in pulling the system toward smaller x menting them in a practical analysis scheme is not trivial, as values (which we here define as the forward direction). Rather, they suggest a segmentation of the trajectory, based on the they should be used in calculating the average reverse work. direction of evolution of the system during each time-step. The rationale for this is that the reversible work portion of the Subsequently, one needs to patch the right segments together work recorded during each red segment equals the negative of to calculate ⟨W⟩F and ⟨W⟩R. Describing one such scheme, the free energy difference for that segment (in the forward together with a method for estimating the associated direction), while its dissipative work is (most likely) still a uncertainties, is what we turn to next. positive value, for the reason discussed above. Adding the work Bin-passing Scheme for Calculation of ⟨W⟩F and ⟨W⟩R. done during the red-segment time-steps in the estimation of the The idea behind the bin-passing scheme for recording ⟨W⟩F forward work, thus cancels some of the reversible work and ⟨W⟩R is that forward (reverse) work should be accounted recorded during the time-steps that actually proceeded in the for when system actually proceeds in the forward (reverse) forward direction (dark blue time-steps), while (most likely) direction at the time-step level. In simulations, when we can 14205 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article check the direction, the system evolves during each time-step, with a value exactly equal to a bin width). Then the temporary there is no reason for mixing together all the back and forth work and path-length variables of the bin are reset to zero and fluctuations that over time progress the system in a certain the algorithm proceeds to reading in the values for the next targeted direction. Rather, we can benefit from our knowledge time-step. of the direction of progress of the system at the smallest (i.e., This scheme involves no greater degree of approximation time-step) level, and then add the work done during each time- than the bin-crossing method, because the accumulation of step to the total amount of work (for the corresponding bin) work within bins in the bin-crossing method already rests upon done in the same direction as the system has progressed during the approximation that (for small enough bins) the force due to that time-step. Simply put, forward work is recorded in a bin surroundings is the same everywhere within the bin. Therefore, only from time-steps during which system actually proceeds each bin passing that occurs fully within a bin measures, to forward, and similarly for the reverse work. within the limits of this approximation, the same reversible Based on this idea, we propose the following scheme for work as the bin crossing method. Data from the time-steps calculating ⟨W⟩F and ⟨W⟩R: perform an SMD FR run, between when a bin-boundary is crossed are discarded, with no the desired initial and final macrostates of the system, noticeable change to the forward or reverse work averages. distinguished by values xA and xB of the reaction coordinate. The extra complexity required in the algorithm for splitting the For each time-step, record (perhaps in a text or binary file) the contributions from such time-steps is not justifiable, given the initial and final values of the reaction coordinate, and also the small amount of extra data that would be gained. external force applied on the system (for a typical two-body When all the raw position and force data is read this way, for scenario, this involves recording initial and final values of two each bin in each direction we have a value for the total work coordinates and two forces, one for each object). We strongly and a counter variable which shows the corresponding number recommend that the analysis be performed postsimulation, as of bin-passings. These two variables are used to find the average this analysis scheme usually takes a few orders of magnitude forward (reverse) work per bin-passing, ⟨W⟩F (⟨W⟩R). These less processing time than the simulation itself, depending on values are then used in eq 2 to calculate ΔG over each bin, and the size of the system and the simulation time invested. There those values are subsequently summed along the reaction path is no extra overhead due to recording the force and position to give the desired PMF, ΔG(x), over the [xA, xB] range. data during the SMD run, as these data have to be recorded for Using this scheme, one generally achieves much better any alternative analysis scheme. Post-simulation analysis also statistics for calculating average external work values, as the preserves raw data, and various analysis algorithms can then be number of bin-passings is usually at least an order of magnitude applied, without the need for repeating the simulation. larger than the corresponding number of bin-crossings, for a After the simulation, the analysis algorithm is given the initial given bin and in a given direction. These better statistics are and final values of the reaction coordinate, and the number of simply the result of exploiting the fluctuating path the system bins (of equal width, although that can be altered, at the price traverses along the reaction path, and assigning each time-step of further complexity in the algorithm). With this information, segment to the proper direction. The work average values an initial mesh of the reaction path can be furnished, with recorded this way also suffer much less from accumulation of spatial boundaries of each bin determined. Conveniently, the dissipative work and simultaneous cancellation of reversible forward direction is then determined as the direction from the work values between different time-step segments that would given initial to the given final value of the reaction coordinate happen in the bin-crossing scheme. xA → xB, with the opposite direction xA ← xB defined as the As we demonstrate below, the bin-passing scheme leads to reverse direction. This choice of direction is arbitrary, since much faster convergence of the RHS of eq 2 to its LHS. reversing it will not affect the shape of the resultant PMF other Nevertheless, we have to consider the possibility that the than, possibly, a constant overall shift. The analysis algorithm statistics gathered using a bin-passing scheme might suffer from then proceeds to read the raw data (recorded with the format internal correlations of data: many instances of bin-passing in a described above), one time-step at a time. For each bin, a path- given direction might result from a single bin-crossing, and the length variable and a temporary work variable is defined in each recorded work values for successive bin-passings might thus be direction. When a time-step datum is read, if the progress of the correlated, i.e., statistically dependent on each other. While system during that step is in the forward (reverse) direction and using such statistics for calculating average work values is safe, its initial and final values of the reaction coordinate fall within we should be cautious in estimating the uncertainty in ⟨W⟩F the ith bin along the reaction path, the temporary work value of that bin is incremented by the value of the forward (reverse) and ⟨W⟩R values (and thus the uncertainty in the PMF) work done during that time-step, and the forward (reverse) calculated this way. As discussed next, we use the block- path-length variable of the bin is incremented with the absolute averaging scheme44 for a more conservative estimation of error value of the displacement along the reaction path during that values. time-step. This process is repeated for subsequent time-steps, Dealing with Errors. Suppose we divide the reaction path until the path-length variable of a bin in a given direction from xA to xB into N bins, each of length δx = |xA − xB|/N. reaches or exceeds the size of the bin (|xB − xA|/(number of During a simulation, we collect ni,F(R) instances of forward- bins)). This marks the completion of one bin-passing. At this (reverse) bin-passing work for the ith bin along the reaction point the temporary work variable of the bin is normalized to path, and we list all those ni instances of bin-passing work in the the value of that bin’s path-length variable, the total forward time-order that they occurred, {W1, W2, ..., Wni}i,F(R). A (reverse) work variable of the bin is incremented with the value correlation length r can effectively be defined for these ni of the normalized forward (reverse) temporary work, and the samplings, such that the mth and the (m + r)th instances of forward (reverse) work-counter variable is incremented by one forward(reverse) bin-passing work along the bin can be (normalization is often required due to the discretization of the considered as statistically independent (one obviously requires motion, which makes it unlikely that any given bin passing ends m + r < ni). Now if we put our ni samplings into ni/r blocks, 14206 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article each of length r, and use these ni/r work averages, rather than choose it to be zero here. Then the value of the free energy the original ni different values to calculate the uncertainty in function at the end of the ith bin is given by ⟨W⟩F(R) across that bin, we can be confident that our ni/r i samplings are statistically independent and thus we are not xB − xA G(xA + i·δx) = ∑ ΔGa , δx = underestimating ⟨δW⟩F(R) for the given bin. This is the idea a=1 N (7) behind block-averaging (Figure 3). Using the ⟨δsW⟩i,F and ⟨δsW⟩i,R values for each bin, individual ΔGi values calculated from eq 2 will have associated uncertainties δs(ΔGi) = ⟨δsW ⟩i2,F + ⟨δsW ⟩i2,R (8) where again the subscript s serves to remind the block-size dependence of our estimation of δ(ΔGi). Using the δs(ΔGi) values found this way along xA → xB, we can estimate the upper bound for the error in our calculated value of the free energy function by the end of the ith bin, δsG(xA + i·δx) Figure 3. Schematic view of block-averaging for work values Wi i corresponding to successive bin-passings along a given bin in a given direction. A block size of s = 3 is used to re-group Wi values into Wj′ δsG(xA + i·δx) = ∑ [δs(ΔGa)]2 values, where, e.g., W1′ = (W1 + W2 + W3)/3. a=1 (9) This equation results in an increasing trend in δsGi values as we The only remaining difficulty is finding the correlation length proceed from xA to xB. We might now recall that the choice of r for the simulation data in hand. That can be simply done by direction xA → xB for accumulating individual ΔGi values to repeating the analysis several times, using different block-sizes. establish the free energy profile G(x) is arbitrary. Subsequently, A block size of s simply means that s successive values of the increasing trend of associated δsGi values in the xA → xB recorded bin-passing work in a given direction and for a given direction also follows that arbitrary choice of direction. If we do bin are averaged together and used as a single value toward the calculation in eq 9 twice, once from xA to xB and the other finding ⟨W⟩F(R) and its associated uncertainty, ⟨δW⟩F(R). By time in the opposite direction, we will obtain two trends of δsG using successively larger values of s, we reduce the chance of that increase in opposite directions. We can convert these two having correlation between each two adjacent block averages of sets of error values into weight factors for their two associated work in a given direction. This generally results in an increase free energy profiles. Recalling that the two free energy profiles in the calculated ⟨δW⟩F(R), using the ni/s block values, rather are exactly the same, except for an overall constant difference, than ni values (corresponding to a block size of 1). The we can shift one to coincide with the other. Corresponding to maximum ⟨δW⟩F(R) will be obtained when s ≈ r. Using larger xA → xB and xA ← xB directions for establishing the G(x) block sizes generally results in a decrease in the calculated profiles, each value of G(x) can now be considered to have two ⟨δW⟩F(R), because when s > r, some uncorrelated values will be different weight factors,45 ws,xA → xB(x) = (δsG(x)xA → xB)−2 and averaged over in a single block, without accounting for their ws,xA ← xB(x) = (δsG(x)xA ← xB)−2. These weight factors can be genuine scatter, which will subsequently result in smaller combined to give a moderated value for δsG(x): calculated ⟨δW⟩F(R). A practical scheme for finding r can not rely solely on 1 δsG(x) = calculating ⟨δW⟩F(R) with different block sizes for one single bin ws , xA → xB(x) + ws , xA ← xB(x) along the reaction path, as each bin generally exhibits a different trend of correlation among successive bin-passing work values. δsG(x)xA → xB . δsG(x)xA ← xB = Rather, a method to establish the error profile along the (δsG(x)xA → xB )2 + (δsG(x)xA ← xB )2 (10) reaction path is needed. The latter should take into account the accumulation of error along the reaction path as we establish δsG(x) values found this way, generally show a trend of starting the PMF using the estimated PMF difference ΔG values for from smaller values at both ends of the reaction path (xA and individual bins. xB) and increasing gradually toward a maximum in the middle. Using a block-size of s, we establish the ⟨W⟩F and ⟨W⟩R This can be understood from eq 8, and the fact that values for each bin along this reaction path. We then use these δsG(x)xA → xB and δsG(x)xA ← xB values show trends of increase values in eq 2 to obtain ΔG across each bin. With ni,F instances of forward and ni,R instances of reverse bin-passing recorded for in the xA → xB and xA ← xB directions, respectively. Equation the ith bin along the reaction path, we will have ni,F/s values of 10 combines these two trends into a more symmetric pattern, forward bin-passing work (with block-size s), and ni,R/s such removing the arbitrariness of unidirectional accumulation of values in the reverse direction. This statistic can be used to error in both δsG(x)xA → xB and δsG(x)xA ← xB profiles. calculate ⟨δsW⟩i,F and ⟨δsW⟩i,R for the ith bin, where the Error profiles established using eqs 8−10, using successively subscript s denotes the general dependence of ⟨δsW⟩i,F(R) on the larger values of block size s, typically have a single maximum in value of the block-size s used. the middle of the reaction path, as we demonstrate in the next We establish ⟨W⟩F and ⟨W⟩R values with their associated section. The value of this maximum increases with s, until s ⟨δsW⟩i,F(R) uncertainties for all the N bins. The initial value of approaches the correlation length r. As explained above, larger free energy function G at xA is arbitrary and does not affect the block sizes (s > r) usually result in a decrease in values of the accumulation of the error, and for the sake of simplicity we error, for averaging over uncorrelated data without accounting 14207 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article for their scatter. We can use this behavior as our guide in Å/ns, equilibrating the system in that configuration for 0.1 ns, estimating r, by establishing error profiles with increasingly moving the Na+ ion back to (2.1 Å, 0, 0) along the x-axis with a larger values of s, until a maximum is reached. Although steering velocity of −7.9 Å/ns, and equilibrating the system in adopting such an error profile as the measure of uncertainty in that configuration for another 0.1 ns. The Cl− ion is our calculated free energy profile might be an overestimate, the constrained at the origin throughout the production run. error-bars obtained this way for the bin-passing method are Each such FR cycle thus takes 2.2 ns, and 10 cycles are often consistent with the bin-to-bin variation in the resultant performed in each of the two production runs, giving a total of free energy profile. We show examples of this procedure and 22 ns of simulation time for each run. The [2.3 Å, 9.5 Å ] range the resultant G and δsG profiles in the next section. of the resultant trajectory has been used in obtaining the PMFs While the convergence rate of the bin-passing algorithm is from both methods. high, there is no guarantee that for a given system and with a The resultant PMFs are shown in Figure 4, together with converged PMF, the error profiles obtained using successively Timko’s PMF for NaCl dissociation obtained using equilibrium larger block sizes would necessarily exhibit a maximum. This classical molecular dynamics with CHARMM27 force-field47 at can be due to lack of enough sampling to encompass several correlation lengths r worth of work values in each direction, while the sampling has been enough already for calculating the underlying PMF. In such cases, we can exploit our knowledge of the physical characteristics of the system. Assume that we know, e.g., that the xA → xB and xA ← xB paths have each been traversed 10 times in a given FR simulation, and that samplings in each two such successive trips can be regarded as uncorrelated. Then for estimating the error in the PMF it will make no sense to use a block-size so large that it gives less than 10 sampled values for any bin along the reaction path, as this would average uncorrelated samplings of work together as one data point. We have used this intuitive rule of thumb for choosing the proper error profiles in the example of peptide- membrane systems below, while the maximum in error-profiles can be established and used for smaller systems when a sufficiently long sampling has been performed (see next section). ■ EXAMPLE 1: PMF OF INTERACTION OF Na+ AND Cl− IONS IN BULK WATER The dissociation PMF for an Na+ and a Cl− ion in water has been studied by many authors, and we use this simple system as our first example to demonstrate the efficiency and accuracy of the bin-passing method. We present here the results of two similar SMD simulations, both performed on the same equilibrated simulation box. For both simulations we have used the TIP3P water model,46 together with the CHARMM force field (version 2747), and both simulations have been performed using NAMD molecular dynamics software,48 at NPT conditions, with atmospheric pressure (P = 1.00 atm) and physiological temperature (T = 310 K). A single Cl− ion is put (strictly fixed) at the origin of the coordinate system, with an Na+ ion at (2.3 Å, 0, 0). A total of Figure 4. (a) Dissociation PMFs for NaCl in water at T = 310 K and P 3916 water molecules are used to build a (nearly) cubical water = 1.00 atm, using TIP3P46 water model and CHARMM2747 force box of size ∼49 Å on each side around these two atoms. The field, obtained from two SMD simulations, and compared with the two ions are held fixed in place using a constraint, while the classical PMF obtained by Timko et al.,49 using the CHARMM27 energy of the whole system is minimized for 1000 steps and force field, at T = 300 K and P = 1 atm. The PMF by Timko et al. has then the system undergoes an equilibration molecular dynamics been calculated over about half the range of reaction coordinate as that run for 5 ns, using a time-step of 2 fs and under NPT studied here using the FR method. We have used a harmonic potential conditions (P = 1.00 atm, T = 310 K). We use the resultant of stiffness k = 500 kcal/(mol·Å2) to steer the Na+ ion with a steering equilibrated box as the initial configuration for both of the speed of 7.9 Å/ns along the trajectory, while the Cl− ion was simulations in this section. constrained at the origin. Error-bars are not shown to allow During the production FR runs, we employ a harmonic distinguishing the PMF curves clearly. (b) Bin-crossing and bin- passing PMFs are given with their corresponding error bars. The error guiding potential of stiffness k = 500 kcal/(mol·Å2) to perform bars for the bin-passing PMF are obtained using a block size of 80, as the steering. Two similar simulations were performed, using explained in Figure 5. (c) The forward portion (moving from larger to identical steering conditions, but with bin-crossing sampling smaller z-values) of the bin-crossing run presented in (a) and (b) is performed in one simulation and bin-passing during the other. used to construct an estimation of the PMF using Jarzynski’s equality Each FR cycle consists of moving the Na+ ion from (2.1 Å, 0, 0) e−βF = ⟨e−βW⟩. The resultant profile has obviously not converged to the to (10 Å, 0, 0) along the x-axis, with a steering velocity of 7.9 underlying PMF, as a comparison with (a) and (b) reveals. 14208 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article 300 K. It is clear in Figure 4 that, for exactly the same amount phosphocholine (DPPC) bilayer at atmospheric pressure and T of simulation time and under the same simulation conditions, = 323 K. the bin-crossing method gives poor results for the dissociation A membrane patch composed of 80 DPPC molecules was PMF of Na+ and Cl− ions in water, while the PMF obtained built, using the CHARMM-GUI50,51 membrane builder. This from the bin-passing algorithm, agrees almost perfectly with patch was then hydrated with a total of 5137 TIP3P46 water Timko’s result. In part (c) of Figure 4 we have shown the PMF molecules on both faces. The approximate dimensions of the obtained using Jarzynski’s equality on the work samplings resultant simulation box were 65 Å × 68 Å × 103 Å along the x, obtained in the bin-crossing simulations. As Jarzynski’s original y, and z directions, respectively, with the normal to the method is unidirectional, only about half of the data from a membrane surface along the z-direction, and its center of mass given series of FR simulations can be used as input to this situated at the origin of the coordinate system. The resultant method, and it is evident that the resultant profile, shown in system was then equilibrated for 60 ns using NAMD48 and the Figure 4c, has not converged to the underlying PMF. We CHARMM3652 force field, with a time-step of 2.00 fs, under should note that using about the same amount of data (∼11 ns constant temperature and pressure (P = 1.00 atm and T = 323 worth of simulation time), the bin-passing method obtains the K), while all three dimensions of the system were allowed to general features of the NaCl dissociation PMF with a much fluctuate in size to let the system settle to an equilibrated better precision (result not shown) than, e.g., the bin-crossing geometry. During the last 5 ns of this equilibration simulation, FR result shown in Figure 4a,b, using 22 ns worth of simulation the area-per-lipid was 60.57 ± 0.84 Å2, exhibiting ∼1.4% time with the conditions explained above. The Jarzynski result fluctuation around the mean value. This value for the average shown in part (c) of Figure 4 thus only serves the purpose of area per lipid headgroup falls within the 57−71.2 Å2 range of exhibiting the relative inefficiency of this unidirectional method experimentally measured values for headgroup area of DPPC at when applied to a small ensemble of trajectories. 50 °C.53 In Figure 5 we establish the error profiles, as discussed in the Nine frames (microstates) were thus obtained, one each 0.5 last section, for the bin-passing PMF presented in Figure 4, ns along this trajectory, which were subsequently used as the initial states for nine parallel sets of FR production runs. The approximate dimensions of these simulation boxes were 60 Å × 67 Å × 105 Å. Each of these boxes was used as the initial configuration for an FR simulation of the maximum duration of 40 ns (totalling to 302.9 ns worth of simulation time). During the FR runs, a single water molecule, which was held fixed at (0, 0, 31 Å) during the equilibrations, was steered to go through the membrane with an average z drift speed of 13.78 Å /ns for a total of 4.5 ns in each direction, and it was held fixed at the two end points (at z = 31 Å and −31 Å) of the trajectory for 0.5 ns between each two successive trips through the membrane. Each FR cycle of the water molecule thus consisted of a 4.5 ns trip Figure 5. PMF error profiles for the FR simulations of NaCl through the membrane from z = 31 Å to −31 Å, 0.5 ns of rest dissociation in water at 310 K, obtained through the bin-passing at z = −31 Å, a 4.5 ns trip back to z = 31 Å, and a final 0.5 ns of algorithm with increasing values of block sizes. A block size of 80 is resting at z = 31 Å, before the next cycle began, totaling to 10 seen to obtain the largest uncertainty values, and the error bars obtained using this value of block size are thus used in Figure 4 on top ns of simulation time. Our 302.9 ns of simulations thus of the bin-passing PMF. encompass over 30 whole FR cycles, performed on nine similar but uncorrelated systems. No restraint was put on either the water molecule or the membrane along the x or y directions, thus, allowing the restrained water molecule to move freely (as using various block sizes from 1 to 160, and it can be seen that a result of thermal fluctuations) in the x and y directions, while these curves reach their maximum at a block size of 80 and its motion along the z axis was restrained through the protocol decrease thereafter for larger block sizes. We have thus used the described above. This was to ensure that the each two error bars obtained with the block size of 80 on top of the bin- successive trips through the width of the membrane were along passing PMF in Figure 4b. These error bars are about half the different trajectories, thus avoiding any bias in sampling the size (or less) of the bin-crossing error bars shown in the same configurations of the phospholipid molecules in the membrane. figure. Comparison of the size of the PMFs and their associated Throughout these FR simulations, the steered water error-bars together demonstrate superior performance of the molecule has been restrained using a spring constant of 0.5 bin-passing method over the traditional bin-crossing method kcal/(mol·Å2), while the center of mass of the DPPC for this system, to the extent that even traces of a third, very membrane patch was restrained to stay at z = 0, using a spring shallow potential well can be distinguished on the bin-passing of stiffness 500 kcal/(mol·Å2). PMF around 7.5 Å. ■ The PMF obtained (using the bin-passing method) from these series of simulations, is shown in Figure 6, together with EXAMPLE 2: PMF FOR PASSAGE OF WATER the (symmetrized) PMFs obtained for the same system by THROUGH A DPPC MEMBRANE various other groups.41,54−57 The uncertainty profiles (using Passage of various molecules through membranes is a classic several different block sizes) for this system are also shown in problem of biological and industrial importance. As the second Figure 7 and exhibit a maximum for a block size of 100. These example of the application of the bin-passing method, we study largest uncertainty estimations are thus used as the error-bars the passage of water through a 1,2-dipalmitoyl-sn-glycero-3- on the bin-passing PMF in Figure 6. The PMFs reported by 14209 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article ⎛ z2 e βΔG(z) ⎞ −1 P = ⎜⎜ ∫z dz⎟⎟ ⎝ 1 D(z) ⎠ (11) where D(z) denotes the local diffusion coefficient. We used the diffusion coefficient profile obtained for this system by Holland et al.41 (at P = 1.00 atm and T = 323 K), together with our PMF for this system, shown in Figure 7, to calculate the right- hand side of eq 11 and found a permeability value of P = (7.92 ± 0.7) × 10−3 cm/s for water−DPPC at the given pressure and temperature. This value is not in perfect agreement with the experimentally measured value of 2.4 × 10−3 cm/s for this Figure 6. PMF for passage of water through a DPPC phospholipid system at T = 319.15 K59 or 3.00 × 10−4 at T = 343.15 K.58 bilayer, as obtained using the bin passing method from FR simulations Nevertheless, this result is in closer agreement with the (at T = 323 K and P = 1 atm), together with various PMFs for the experimental values than those obtained through previous same system obtained by other groups.41,54−57 Only, the bin-passing simulation studies, as shown in Figure 8. PMF as given here is not symmetrized, and it exhibits a maximum asymmetry of only ∼0.17 kcal/mol. The value z = 0 corresponds to the center of the membrane. Figure 8. Comparison of water permeability of the DPPC bilayer, as Figure 7. Error profiles for the 302.9 ns long FR simulation of the obtained from various studies. The reported values given by Jansen passage of water through a DPPC bilayer. A block size of 100 is seen to and Blume (343.15 K)58 and Lawaczeck (319.15 K)59 are obtained obtain the largest estimated uncertainties in the PMF, and successively experimentally, while the other four results41,55−57 and the value higher block sizes result in reduced estimations of the uncertainty. The obtained in this work are from simulation studies. uncertainty values obtained with a block size of 100 are used on the bin-passing PMF shown in Figure 6. As can be seen from eq 11, permeability depends exponentially on the PMF and only linearly on the diffusion coefficient. The better agreement with the experimental results other groups are obtained using various methods, force-fields, of the permeability obtained using the PMF found through the and simulation times. The longest simulation among these has bin-passing method (compared with those reported from other been the 30 ns run performed by Holland et al.,41 and the simulation studies) thus supports the validity of the larger barrier height obtained here, shown in Figure 6. ■ shortest has perhaps been the classic and pioneering studies by Marrink and Berendsen,55 where they could afford only to spend very few nanoseconds of simulation time (total EXAMPLE 3: PMF OF INTERACTION OF THE simulation time unclear). With the exception of the PMF HHC-36 PEPTIDE WITH TWO MODEL MEMBRANES reported by Shinoda et al.,56 all other groups calculated the As our third example of the application of the bin-passing PMF for only one leaflet and mirrored the result to produce the method, we establish the adsorption PMF for the antimicrobial full profile. Nevertheless, Shinoda et al.56 and Sugii et al.57 did peptide HHC-36 (WRWWKWWRR)43 to two different model not provide error bars on their PMFs, and the degree of membranes. HHC-36 has already been shown to perform as a convergence of their results is thus undetermined. strong antimicrobial peptide, through a series of in vivo and in In Figure 6, the bin-passing PMF shows a deviation from vitro experiments. Specifically, this peptide has been tested symmetry of no more than ∼0.17 kcal/mol, over the whole against strains of multidrug-resistant bacteria, such as reaction path of length 60 Å. The PMFs obtained by the other Pseudomonas aeruginosa and Staphylococcus aureus and has groups were symmetrized (by amounts not reported), due to surpassed regularly used antibiotics in antimicrobial activity. incomplete sampling. Our bin-passing PMF for this system HHC-36 has also shown minimal toxicity toward host exhibits a larger barrier height (by ∼1−2 kcal/mol) compared mammalian cells, specifically toward red blood cells.43 While to those obtained by other groups. However, our PMF results these experiments can not indicate the precise mechanism of in a better agreement with experimentally measured perme- interaction of this peptide against either mammalian or bacterial ability values when used to obtain the DPPC bilayer cell membranes, they do suggest that the peptide adsorbs less permeability to water at the given temperature, using the strongly to mammalian cell membranes and more so to the equation55 bacterial cell membranes. This suggestion is supported by the 14210 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article fact that mammalian cell membranes are rich in zwitterionic above the membrane surface, with the z-component of its phospholipids, while the cytoplasmic membrane of bacteria is center of mass at 55 Å, that is, 55 Å along the direction normal rich in anionic phospholipids. With five units of positive electric to the membrane surface. The peptide was then moved toward charge, one would expect this peptide to show more attraction the membrane with a drift velocity of 5 Å/ns, along the z-axis, to negatively charged phospholipids like phosphatidylglycerol for 7 ns, until the z-coordinate of its center of mass reached 20 (PG), compared to the zwitterionic phosphatidylcholine (PC). Å. We kept the peptide there for half a nanosecond to let the Strains of methicillin-resistant Staphylococcus aureus (a Gram- system settle down, and then steered the system along the positive bacterium) contain between 70 to over 85% of the PG opposite path, moving peptide back along the z-axis to z = 55 in their total phospholipid content.60 PG also comprises 18% of Å, followed by another 0.5 ns of relaxation, which marked the the phospholipid content of the Salmonella typhimurium, and completion of one round of an FR run with |vF| = |vR| = 5 Å/ns. 19% of that of Escherichia coli (both Gram-negative), with the A snapshot of the initial configuration of this cycle, together major phospholipid component in both cases being phospha- with a depiction of distances described above, is shown in tidylethanolamine (PE), which comprises about 75% of the Figure 9, where the water molecules are shown very faintly by phospholipid content of the membrane of Salmonella dotted lines to help keep the peptide and the membrane surface typhimurium, and 69% in the case of Escherichia coli.61 PC is visible in the figure. the major phospholipid component of mammalian erythro- cytes. While it comprises about 30% of the phospholipid content of the human erythrocyte membrane, 62,63 its distribution is not symmetric between inner and outer leaflets: the outer leaflet is has about 45% PC content among its phospholipids, while the inner leaflet has only 14%.62 To study the interaction of HHC-36 peptide with mammalian and bacterial membranes we prepared two model membrane patches, one composed of 128 molecules of the 1- palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) mole- cules, with 64 phospholipid molecules per leaflet, and the other with a 3:1 mixture of 1-palmitoyl-2-oleoyl-sn-glycero-3- phosphoethanolamine (POPE) to 1-palmitoyl-2-oleoyl-sn-glyc- ero-3-phosphoglycerol (POPG) molecules, using 32 POPG and 96 POPE molecules, with 48 POPE and 16 POPG molecules in each leaflet. Both membranes were built using the CHARMM- Figure 9. General scheme for the SMD simulations of the HHC-36 adsorption to the POPC and POPE/POPG membrane patches. GUI membrane builder50,51 and then hydrated in VMD Center of mass of the peptide is moved along the z-axis from z = 55 to molecular visualization software.64 The resultant simulation 20 Å and back, with an average steering speed of 5.00 Å/ns, while the boxes were then merged each with a separate box containing a center of mass of the membrane patch is restrained at z = 0. hydrated HHC-36 molecule, obtained from over 600 ns refolding simulation, initiated from a refolded structure of the peptide molecule. All of these simulations (initial equilibration, For steering, the HHC-36 molecule against the POPE/ refolding simulation of the stretched peptide structure, final POPG membrane patch in the fashion just described, we use a equilibration, and the production runs) have been done using harmonic restraint of kp = 500 kcal/(mol·Å2), with a harmonic NAMD48 at T = 310 K and P = 1.00 atm, using TIP3P water restraint of km = 10000 kcal/(mol·Å2) holding the z-coordinate model,46 with CHARMM2747 force field for water molecules, of the center of mass of the membrane essentially fixed at zero. ions, and the protein, and CHARMM3652 force field for the For the peptide-POPC system, these values were kp = 100 kcal/ membrane patches. Proper numbers of sodium and chloride (mol·Å2) and km = 5000 kcal/(mol·Å2), respectively. A total of ions have been added to each simulation box to obtain near 743.87 ns of production FR run was collected for the peptide- physiological NaCl concentration (∼0.15 M), and also to POPE/POPG system, and 611 ns for the peptide-POPC render the simulation boxes electrically neutral. The final system. The collected force and positions data from all the POPE/POPG box (containing the peptide) has been time-steps for each system were analyzed by the bin-passing equilibrated for 71.7 ns, and the POPC box (also containing algorithm described in the earlier sections, and the resultant the peptide) for 175.6 ns. POPC by itself does not form very adsorption PMFs are shown in Figure 10. One can see in that stable membranes, and thus, we performed a longer figure that the PMF well for the peptide’s interaction with the equilibration run for the POPC box, compared to the one charged membrane (POPE−POPG) is much deeper compared containing POPE/POPG. to that of the zwitterionic membranes (POPC). The well- From the last 20 ns of the equilibration run of the POPE/ depths are 13.41 and 4.75 kcal/mol for the POPE−POPG and POPG box, 11 frames have been picked, at 2 ns intervals, which POPC membranes, respectively. Both curves are essentially flat were then used as the starting configurations of the 11 different in the 43−47 Å range, which means that the peptide is out of FR SMD peptide-membrane runs. Similarly, from the last 24 ns the effective interaction range of the membranes at this distance of the POPC+peptide box, 13 frames have been extracted at interval, and we have shifted both curves to have an average of intervals 2 ns apart, and these have been used as the initial zero in this region. We wish to emphasize the efficiency of the configurations for 13 different FR SMD runs. bin-passing method in obtaining these curves, giving Each membrane has its interfaces parallel to the x−y plane, satisfactory PMFs within less than a microsecond of simulation and the z coordinate of its center of mass has been kept close to time. Applying the bin-crossing method would require much zero by means of harmonic restraints. In our FR simulations longer simulation times to achieve a converged PMF for either with both types of membranes, the peptide was initially held of the two peptide-membrane systems presented here. 14211 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article conveniently chosen fast density functional theory (DFT) methods in between the pure classical steps, which are in turn influenced by the outcome of quantum calculations. This way a practical distinction is made between the classical and quantum time-scales (following the same idea presented in the Born− Oppenheimer approximation67), which makes it possible to apply non-equilibrium work theorems to the classical steps. Due to the fact that the total distance is accumulated to determine the number of bin-passings (as opposed to the total displacement used in the bin-crossing method), a much higher number of work samples is obtained in this method, further improving the statistics. The interrelation among those contributions can be vividly seen in the central formulae of the FR method, and can be effectively implemented when the Figure 10. PMFs for adsorption of HHC-36 molecule to a POPC (pink curve) and a POPE/POPG (dark blue curve) membrane patch, notion of dissipative work is considered within the framework both at T = 310 K and P = 1.00 atm. TIP3P water model has been of the second law of thermodynamics. These considerations are used, together with the CHARMM2747 force field for the peptide, fundamental to the design of the bin-passing method presented water, and ions, and the CHARMM3652 force field for the membranes. here, and contribute to its higher efficiency compared to the A drift velocity of 5 Å/ns has been used for steering the peptide more elementary bin-crossing scheme. Our experimentation toward each membrane surface and away from it. Error-bars are with the bin-crossing method (results shown for the NaCl obtained using block-sizes of 400 for the POPE/POPG system and dissociation in water) makes it clear that this method performs 300 for the POPC system, as these where the largest block-size values poorly in obtaining the PMFs for thermodynamic systems. This that for each bin gave at least as many samplings as the bin-crossing is to be expected from the arguments presented above: scheme would obtain (∼50 for the POPE/POPG system and ∼41 for accumulation of the dissipative work and cancellation of the POPC system.). Center of mass of each membrane patch is restrained at z = 0. reversible work contributions in both the forward and the reverse work averages, when sampling is done through the bin- crossing method, results in poor convergence to the PMF, and Furthermore, as we discussed earlier, the inherent deficiencies larger uncertainties in the resultant profiles. In contrast, the bin- of the bin-crossing method in obtaining reliable average passing method exploits the thermal fluctuations occurring forward and reverse work values might make it impossible to during SMD simulations, and results in highly improved obtain converged PMFs using that method, even from multi- convergence and smaller error bars. microsecond simulations of these relatively large systems. The convergence rate generally decreases with the size of the The larger well-depth of the PMF for the peptide’s system under study, as is the case with any other PMF interaction with the POPE−POPG membrane, compared to calculation method. Nevertheless, we showed that a sub- that with the POPC membrane, supports the higher affinity of microsecond long simulation is sufficient for obtaining the peptide to adsorb to the POPE−POPG membrane. satisfactory PMFs for systems as large as a hydrated peptide- Recalling that POPE−POPG is a memetic of bacterial membrane pair. Furthermore, the block-size variation scheme cytoplasmic membrane and POPC represents a model for that we present and discuss puts a tight bound on the maximum mammalian cell membranes, this observation is in agreement uncertainty of the PMFs obtained using the bin-passing with in vivo and in vitro experiments43 and also isothermal method. titration calorimetry experiments65 that have collectively We suspect that lack of a proper analysis scheme and established HHC-36 as a strong antimicrobial agent with software has been a major reason for the FR method33 not minimal hemolytic activity. becoming as widely employed by simulators as (we believe) it ■ CONCLUSION Based on the FR method,33,38,39 the bin-passing scheme was should be. The analysis scheme presented here can be performed using the bin-passing analyzer software that we have made available as a free software at https://github.com/ shown here to effectively and precisely obtain the PMFs for 1particle/bin-passing_analyzer under the terms of the GNU systems of varying degrees of complexity. The key point in General Public License (version 3). This software has been devising this scheme and thus to obtain such free energy used for performing the post simulation analyses for all the profiles is to carefully consider the reversible and dissipative PMFs presented in this paper. Considering the remarkable work contributions during each time-step of the steered convergence rate and precision of the bin-passing method and molecular dynamics simulations, and to assign them accord- the advantages it provides (such as manual post-simulation ingly to forward and reverse directions. This is in contrast to binning of the reaction path), we hope that these methods will the bin-crossing scheme that relies on cancellation of be adopted by the community of chemical and biomolecular simulators. ■ undesirable work contributions merely through the averaging process, but accumulates the error of both the forward and reverse dissipative work in each unidirectional bin. This method AUTHOR INFORMATION can be generalized to mixed quantum mechanical-molecular Corresponding Authors mechanical (QM/MM) simulations, such as those invoking the *E-mail:

. Car−Parrinello scheme.66 some groups have demonstrated the *E-mail:

. application of Jarzynski’s equality to calculate PMFs from *E-mail:

steered QM/MM simulations.30,31 The quantum mechanical Notes variations of the electronic structure are often handled by The authors declare no competing financial interest. 14212 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B ■ Article ACKNOWLEDGMENTS (19) Cuendet, M. A.; Michielin, O. Protein−protein interaction investigated by steered molecular dynamics: the TCR−pMHC The authors are grateful to the Natural Sciences and complex. Biophys. J. 2008, 95, 3575−3590. Engineering Research Council of Canada for financial support. (20) Boechi, L.; Martí, M. A.; Vergara, A.; Sica, F.; Mazzarella, L.; Computational resources for the simulations presented in this Estrin, D. A.; Merlino, A. Protonation of histidine 55 affects the work were provided by SHARCNET, SciNet,68 and Calcul oxygen access to heme in the alpha chain of the hemoglobin from the Québec HPC consortia of Compute Canada. Antarctic fish Trematomus bernacchii. IUBMB Life 2011, 63, 175−182. ■ REFERENCES (1) Jarzynski, C. Nonequilibrium equality for free energy differences. (21) Morzan, U. N.; Capece, L.; Marti, M. A.; Estrin, D. A. Quaternary structure effects on the hexacoordination equilibrium in rice hemoglobin rHb1: Insights from molecular dynamics simulations. Phys. Rev. Lett. 1997, 78, 2690−2693. Proteins 2013, 81, 863−873. (2) Jarzynski, C. Equilibrium free-energy differences from non- (22) Patel, J. S.; Berteotti, A.; Ronsisvalle, S.; Rocchia, W.; Cavalli, A. equilibrium measurements: A master-equation approach. Phys. Rev. E Steered molecular dynamics simulations for studying protein-ligand 1997, 56, 5018−5035. interaction in cyclin-dependent kinase 5. J. Chem. Inf. Model. 2014, 54, (3) Crooks, G. Path-ensemble averages in systems driven far from 470−480. equilibrium. Phys. Rev. E 2000, 61, 2361−2366. (23) Yu, T.; Lee, O. S.; Schatz, G. C. Steered molecular dynamics (4) Hummer, G.; Szabo, A. Free energy reconstruction from studies of the potential of mean force for peptide amphiphile self- nonequilibrium single-molecule pulling experiments. Proc. Natl. Acad. assembly into cylindrical nanofibers. J. Phys. Chem. A 2013, 117, Sci. U.S.A. 2001, 98, 3658−3661. 7453−7460. (5) Liphardt, J.; Dumont, S.; Smith, S. B.; Tinoco, I.; Bustamante, C. (24) Hwang, H.; Schatz, G. C.; Ratner, M. A. Steered molecular Equilibrium information from nonequilibrium measurements in an dynamics studies of the potential of mean force of a Na+ or K+ ion in experimental test of Jarzynski’s equality. Science 2002, 296, 1832− a cyclic peptide nanotube. J. Phys. Chem. B 2006, 110, 26448−26460. 1835. (25) L. Boechi, L.; Martí, M. A.; Milani, M.; Bolognesi, M.; Luque, F. (6) Yangyuoru, P. M.; Dhakal, S.; Yu, Z.; Koirala, D.; Mwongela, S. J.; Estrin, D. A. Structural determinants of ligand migration in M.; Mao, H. Single-molecule measurements of the binding between Mycobacterium tuberculosis truncated hemoglobin O. Proteins 2008, small molecules and DNA aptamers. Anal. Chem. 2012, 84, 5298− 73, 372−379. 5303. (26) Caballero, J.; Zamora, C.; Aguayo, D.; Yanez, C.; González-Nilo, (7) Park, S.; Schulten, K. Calculating potentials of mean force from F. D. Study of the interaction between progesterone and β- steered molecular dynamics simulations. J. Chem. Phys. 2004, 120, cyclodextrin by electrochemical techniques and steered molecular 5946−5961. dynamics. J. Phys. Chem. B 2008, 112, 10194−10201. (8) Minh, D.; McCammon, J. Springs and speeds in free energy (27) Jensen, M. Ø.; Yin, Y.; Tajkhorshid, E.; Schulten, K. Sugar reconstruction from irreversible single-molecule pulling experiments. J. transport across lactose permease probed by steered molecular Phys. Chem. B 2008, 112, 5892−5897. dynamics. Biophys. J. 2007, 93, 92−102. (9) Jensen, M.; Park, S.; Tajkhorshid, E.; Schulten, K. Energetics of (28) Boechi, L.; Mañez, P. A.; Luque, F. J.; Marti, M. A.; Estrin, D. A. glycerol conduction through aquaglyceroporin GlpF. Proc. Natl. Acad. Unraveling the molecular basis for ligand binding in truncated Sci. U.S.A. 2002, 99, 6731−6736. hemoglobins: The trHbO Bacillus subtilis case. Proteins 2010, 78, 962− (10) Xiong, H.; Crespo, A.; Marti, M.; Estrin, D.; Roitberg, A. Free 970. energy calculations with non-equilibrium methods: applications of the (29) Wang, Y.; Schulten, K.; Tajkhorshid, E. What makes an Jarzynski relationship. Theor. Chem. Acc. 2006, 116, 338−346. aquaporin a glycerol channel? A comparative study of AqpZ and GlpF. (11) Forti, F.; Boechi, L.; Estrin, D. A.; Marti, M. A. Comparing and Structure 2005, 13, 1107−1118. combining implicit ligand sampling with multiple steered molecular (30) Raugei, S.; Cascella, M.; Carloni, P. A proficient enzyme: dynamics to study ligand migration processes in heme proteins. J. insights on the mechanism of orotidine monophosphate decarboxylase Comput. Chem. 2011, 32, 2219−2231. from computer simulations. J. Am. Chem. Soc. 2004, 126, 15730− (12) Jensen, M. Ø.; Park, S.; Tajkhorshid, E.; Schulten, K. Energetics 15737. of glycerol conduction through aquaglyceroporin GlpF. Proc. Natl. (31) Crespo, A.; Martí, M. A.; Estrin, D. A.; Roitberg, A. E. Multiple- Acad. Sci. U.S.A. 2002, 99, 6731−6736. steering QM-MM calculation of the free energy profile in chorismate (13) Gupta, A. N.; Vincent, A.; Neupane, K.; Yu, H.; Wang, F.; mutase. J. Am. Chem. Soc. 2005, 127, 6940−6941. Woodside, M. T. Experimental validation of free-energy-landscape (32) Izrailev, S.; Stepaniants, S.; Isralewitz, B.; Kosztin, D.; Lu, H.; reconstruction from non-equilibrium single-molecule force spectros- Molnar, F.; Wriggers, W.; Schulten, K. Lecture Notes in Computa- copy measurements. Nat. Phys. 2011, 7, 631−634. (14) Moreno, D. M.; Martí, M. A.; Biase, P. M. D.; Estrin, D. A.; tional Science and Engineering. In Computational Molecular Dynamics: Demicheli, V.; Radi, R.; Boechi, L. Exploring the molecular basis of Challenges, Methods, Ideas - Proceedings of the 2nd International human manganese superoxide dismutase inactivation mediated by Symposium on Algorithms for Macromolecular Modelling; Deuflhard, tyrosine 34 nitration. Arch. Biochem. Biophys. 2011, 507, 304−309. P., Hermans, J., Leimkuhler, B., Mark, A., Reich, S., Skeel, R., Eds.; (15) Bidon-Chanal, A.; Martí, M. A.; Crespo, A.; Milani, M.; Orozco, Springer-Verlag: Berlin, 1998; Vol. 4; Chapter “Steered Molecular M.; Bolognesi, M.; Luque, F. J.; Estrin, D. A. Ligand-induced Dynamics”, pp 39−65. dynamical regulation of NO conversion in Mycobacterium tuber- (33) Kosztin, I.; Barz, B.; Janosi, L. Calculating potentials of mean culosis truncated hemoglobin-N. Proteins: Struct., Funct., Bioinf. 2006, force and diffusion coefficients from nonequilibrium processes without 64, 457−464. Jarzynski’s equality. J. Chem. Phys. 2006, 124, 064106. (16) Kang, Y.; Liu, Y. C.; Wang, Q.; Shen, J. W.; Wu, T.; Guan, W. J. (34) Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free On the spontaneous encapsulation of proteins in carbon nanotubes. energy calculation from steered molecular dynamics simulations using Biomaterials 2009, 30, 2807−2815. Jarzynski’s equality. J. Chem. Phys. 2003, 119, 3559−3566. (17) Zhang, D.; Gullingsrud, J.; McCammon, J. A. Potentials of mean (35) Jarzynski, C. Rare events and the convergence of exponentially force for acetylcholine unbinding from the alpha7 nicotinic acetylcho- averaged work values. Phys. Rev. E 2006, 73, 046105. line receptor ligand-binding domain. J. Am. Chem. Soc. 2006, 128, (36) Collin, D.; Ritort, F.; Jarzynski, C.; Smith, S. B.; T, I., Jr; 3019−3026. Bustamante, C. Verification of the Crooks fluctuation theorem and (18) Martí, M. A.; Lebrero, M. C. G.; Roitberg, A. E.; Estrin, D. A. recovery of RNA folding free energies. Nature 2005, 437, 231−234. Bond or cage effect: how nitrophorins transport and release nitric (37) Pohorille, A.; Jarzynski, C.; Chipot, C. Good practices in free- oxide. J. Am. Chem. Soc. 2008, 130, 1611−1618. energy calculations. J. Phys. Chem. B 2010, 114, 10235−10253. 14213 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214 The Journal of Physical Chemistry B Article (38) Forney, M.; Janosi, L.; Kosztin, I. Calculating free-energy (59) Lawaczeck, R. On the permeability of water molecules across profiles in biomolecular systems from fast nonequilibrium processes. vesicular lipid bilayers. J. Membr. Biol. 1979, 51, 229−261. Phys. Rev. E 2008, 78, 051913. (60) Mishra, N. N.; Yang, S. J.; Sawa, A.; Rubio, A.; Nast, C. C.; (39) Fabritiis, G. D.; Coveney, P.; Villa‘-Freixa, J. Energetics of K+ Yeaman, M. R.; Bayer, A. S. Analysis of cell membrane characteristics permeability through Gramicidin A by forward-reverse steered of in vitro-selected daptomycin-resistant strains of methicillin-resistant molecular dynamics. Proteins 2008, 73, 185−194. Staphylococcus aureus. Antimicrob. Agents Chemother. 2009, 53, 2312− (40) NategholEslam, M.; Holland, B. W.; Gray, C. G.; Tomberli, B. 2318. Drift-oscillatory steering with the forward-reverse method for (61) Ames, G. F. Lipids of Salmonella typhimurium and Escherichia calculating the potential of mean force. Phys. Rev. E 2011, 83, 021114. coli: structure and metabolism. J. Bacteriol. 1968, 95, 833−843. (41) Holland, B. W.; Gray, C. G.; Tomberli, B. Calculating diffusion (62) Dodge, J. T.; Phillips, G. B. Composition of phospholipids and and permeability coefficients with the oscillating forward-reverse of phospholipid fatty acids and aldehydes in human red cells. J. Lipid method. Phys. Rev. E 2012, 86, 036707. Res. 1967, 8, 667−675. (42) Jarzynski, C. Progress in Mathematical Physics; IEqualities and (63) Virtanen, J. A.; Cheng, K. H.; Somerharju, P. Phospholipid Inequalities: Irreversibility and the Second Law of Thermodynamics at composition of the mammalian red cell membrane can be rationalized the Nanoscale. In Time; Duplantier, B., Ed.; 2013; Vol. 63, pp 145− by a superlattice model. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 4964− 172, 4969. (43) Cherkasov, A.; Hilpert, K.; Jenssen, H.; Fjell, C. D.; Waldbrook, (64) Humphrey, W.; Dalke, A.; Schulten, K. VMD - visual molecular M.; Mullaly, S. C.; Volkmer, R.; Hancock, R. E. W. Use of artificial dynamics. J. Mol. Graph. 1996, 14, 33−38. intelligence in the design of small peptide antibiotics effective against a (65) Nichols, M.; Kuljanin, M.; Nategholeslam, M.; Hoang, T.; broad spectrum of highly antibiotic-resistant superbugs. ACS Chem. Vafaei, S.; Tomberli, B.; Gray, C. G.; DeBruin, L.; Jelokhani-Niaraki, Biol. 2008, 4, 65−74. M. Dynamic turn conformation of a short tryptophan-rich cationic (44) Rapaport, D. C. The Art of Molecular Dynamics Simulation, 2nd antimicrobial peptide and its interaction with phospholipid mem- ed.; Cambridge University Press: New York, 2004; p 86. branes. J. Phys. Chem. B 2013, 117, 14697−14708. (45) Taylor, J. R. An Introduction to Error Analysis: the Study of (66) Car, R.; Parrinello, M. Unified approach for molecular dynamics Uncertainties in Physical Measurements; University Science Books: and density-functional theory. Phys. Rev. Lett. 1985, 55, 2471−2474. Herndon, VA, 1997; Chapter 7. (67) Born, M.; Oppenheimer, R. Zur quantentheorie der molekeln. (46) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. Ann. Phys. 1927, 389, 457−484. W.; Klein, M. L. Comparison of simple potential functions for (68) Loken, C.; et al. SciNet: lessons learned from building a power- simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. efficient top-20 system and data centre. J. Phys.: Conf. Ser. 2010, (47) Foloppe, N.; MacKerell, A. D., Jr All-atom empirical force field 012026. for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data. J. Comput. Chem. 2000, 21, 86−104. (48) Phillips, J. .; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable molecular dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (49) Timko, J.; Bucher, D.; Kuyucak, S. Dissociation of NaCl in water from ab initio molecular dynamics simulations. J. Chem. Phys. 2010, 132, 114510. (50) Jo, S.; Kim, T.; Im, W. Automated builder and database of protein/membrane complexes for molecular dynamics simulations. PLoS One 2007, 2, e880. (51) Jo, S.; Lim, J. B.; Klauda, J. B.; Im, W. CHARMM-GUI membrane builder for mixed bilayers and its application to yeast membranes. Biophys. J. 2009, 97, 50−58. (52) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D., Jr CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J. Comput. Chem. 2010, 31, 671−690. (53) Nagle, J. F.; Tristram-Nagle, S. Structure of lipid bilayers. Biochim. Biophys. Acta, Biomembr. 2000, 1469, 159−195. (54) Bemporad, D.; Essex, J. W.; Luttmann, C. Permeation of small molecules through a lipid bilayer: a computer simulation study. J. Phys. Chem. B 2004, 108, 4875−4884. (55) Marrink, S. J.; Berendsen, H. J. C. Simulation of water transport through a lipid membrane. J. Phys. Chem. 1994, 98, 4155−4168. (56) Shinoda, W.; Mikami, M.; Baba, T.; Hato, M. Molecular dynamics study on the effects of chain branching on the physical properties of lipid bilayers: 2. Permeability. J. Phys. Chem. B 2004, 108, 9346−9356. (57) Sugii, T.; Takagi, S.; Matsumoto, Y. A molecular-dynamics study of lipid bilayers: effects of the hydrocarbon chain length on permeability. J. Chem. Phys. 2005, 123, 184714. (58) Jansen, M.; Blume, A. A comparative study of diffusive and osmotic water permeation across bilayers composed of phospholipids with different head groups and fatty acyl chains. Biophys. J. 1995, 68, 997−1008. 14214 dx.doi.org/10.1021/jp504942t | J. Phys. Chem. B 2014, 118, 14203−14214