Single-ensemble nonequilibrium path-sampling estimates of free energy differences

Nonuniform convergence in moment expansions of integral work relations

Exponential averages that appear in integral fluctuation theorems can be recast as a sum over moments of thermodynamic observables. We use two examples to show that such moment series can exhibit nonuniform convergence in certain singular limits. The first example is a simple model of a process with measurement and feedback. In this example, the limit of interest is that of error-free measurements. The second system we study is an ideal gas particle inside an (infinitely) fast expanding piston. Both examples show qualitative similarities; the low-order moments are close to their limiting value, while high-order moments strongly deviate from their limit. As the limit is approached the transition between the two groups of moments is pushed toward higher and higher moments. Our findings highlight the importance of the ordering of limits in certain nonequilibrium-related calculations.

Enhanced Jarzynski free energy calculations using weighted ensemble

The free energy of transitions between stable states is the key thermodynamic quantity that governs the relative probabilities of the forward and reverse reactions and the ratio of state probabilities at equilibrium. The binding free energy of a drug and its receptor is of particular interest, as it serves as an optimization function for drug design. Over the years, many computational methods have been developed to calculate binding free energies, and while many of these methods have a long history, issues such as convergence of free energy estimates and the projection of a binding process onto order parameters remain. Over 20 years ago, the Jarzynski equality was derived with the promise to calculate equilibrium free energies by measuring the work applied to short nonequilibrium trajectories. However, these calculations were found to be dominated by trajectories with low applied work that occur with extremely low probability. Here, we examine the combination of weighted ensemble algorithms with the Jarzynski equality. In this combined method, an ensemble of nonequilibrium trajectories are run in parallel, and cloning and merging operations are used to preferentially sample low-work trajectories that dominate the free energy calculations. Two additional methods are also examined: (i) a novel weighted ensemble resampler that samples trajectories directly according to their importance to the work of work and (ii) the diffusion Monte Carlo method using the applied work as the selection potential. We thoroughly examine both the accuracy and efficiency of unbinding free energy calculations for a series of model Lennard-Jones atom pairs with interaction strengths ranging from 2 kcal/mol to 20 kcal/mol. We find that weighted ensemble calculations can more efficiently determine accurate binding free energies, especially for deeper Lennard-Jones well depths.

Best Practices for Alchemical Free Energy Calculations [Article v1.0]

Alchemical free energy calculations are a useful tool for predicting free energy differences associated with the transfer of molecules from one environment to another. The hallmark of these methods is the use of "bridging" potential energy functions representing alchemical intermediate states that cannot exist as real chemical species. The data collected from these bridging alchemical thermodynamic states allows the efficient computation of transfer free energies (or differences in transfer free energies) with orders of magnitude less simulation time than simulating the transfer process directly. While these methods are highly flexible, care must be taken in avoiding common pitfalls to ensure that computed free energy differences can be robust and reproducible for the chosen force field, and that appropriate corrections are included to permit direct comparison with experimental data. In this paper, we review current best practices for several popular application domains of alchemical free energy calculations performed with equilibrium simulations, in particular relative and absolute small molecule binding free energy calculations to biomolecular targets.

Transient probability currents provide upper and lower bounds on non-equilibrium steady-state currents in the Smoluchowski picture

Probability currents are fundamental in characterizing the kinetics of nonequilibrium processes. Notably, the steady-state current Jss for a source-sink system can provide the exact mean-first-passage time (MFPT) for the transition from the source to sink. Because transient nonequilibrium behavior is quantified in some modern path sampling approaches, such as the "weighted ensemble" strategy, there is strong motivation to determine bounds on Jss-and hence on the MFPT-as the system evolves in time. Here, we show that Jss is bounded from above and below by the maximum and minimum, respectively, of the current as a function of the spatial coordinate at any time t for one-dimensional systems undergoing overdamped Langevin (i.e., Smoluchowski) dynamics and for higher-dimensional Smoluchowski systems satisfying certain assumptions when projected onto a single dimension. These bounds become tighter with time, making them of potential practical utility in a scheme for estimating Jss and the long time scale kinetics of complex systems. Conceptually, the bounds result from the fact that extrema of the transient currents relax toward the steady-state current.

Trajectory stratification of stochastic dynamics

We present a general mathematical framework for trajectory stratification for simulating rare events. Trajectory stratification involves decomposing trajectories of the underlying process into fragments limited to restricted regions of state space (strata), computing averages over the distributions of the trajectory fragments within the strata with minimal communication between them, and combining those averages with appropriate weights to yield averages with respect to the original underlying process. Our framework reveals the full generality and flexibility of trajectory stratification, and it illuminates a common mathematical structure shared by existing algorithms for sampling rare events. We demonstrate the power of the framework by defining strata in terms of both points in time and path-dependent variables for efficiently estimating averages that were not previously tractable.

Calculating the free energy difference by applying the Jarzynski equality to a virtual integrable system

The Jarzynski equality (JE) provides a nonequilibrium method to measure and calculate the free energy difference (FED). Note that if two systems share the same Hamiltonian at two equilibrium states, respectively, they share the same FED between these two equilibrium states as well. Therefore the calculation of the FED of a system may be facilitated by considering instead another virtual system designed to this end. Taking advantage of this flexibility and the JE, we show that by introducing an integrable virtual system, the evolution problem involved in the JE can be solved. As a consequence, FED is expressed in the form of an equilibrium equality, in contrast with the nonequilibrium JE it is based on. Numerically, this result allows FED to be computed by sampling the canonical ensemble directly and the computational cost can be significantly reduced. The effectiveness and efficiency of this scheme are illustrated with numerical studies of several representative model systems.

Jarzynski matrix equality: Calculating the free-energy difference by nonequilibrium simulations with an arbitrary initial distribution

The Jarzynski equality (JE) method, which relates the work of a nonequilibrium process to the free-energy difference between its initial and final states, provides an efficient way to calculate free energies of thermodynamic systems in simulations or experiments. However, more extensive applications of the JE are hindered by the requirement that the initial state must be in equilibrium. In this work we extend the JE method to be the Jarzynski matrix equality (JME) method, which relates the work of trajectories connecting metastable conformational regions to their local free energies, and thus we can estimate the free energy from the nonequilibrium trajectories starting from an almost arbitrary initial distribution. We then apply the JME to toy models, Lennard-Jones fluids, and polymer chain models, demonstrating its efficiency in free-energy calculations with satisfactory accuracy. The JME extends the applicability of the nonequilibrium methods to complex systems whose initial equilibrium states are difficult to reach.

Elastic Barrier Dynamical Freezing in Free Energy Calculations: A Way To Speed Up Nonequilibrium Molecular Dynamics Simulations by Orders of Magnitude

An important issue concerning computer simulations addressed to free energy estimates via nonequilibrium work theorems, such as the Jarzynski equality [Phys. Rev. Lett. 1997, 78, 2690], is the computational effort required to achieve results with acceptable accuracy. In this respect, the dynamical freezing approach [Phys. Rev. E 2009, 80, 041124] has been shown to improve the efficiency of this kind of simulations, by blocking the dynamics of particles located outside an established mobility region. In this report, we show that dynamical freezing produces a systematic spurious decrease of the particle density inside the mobility region. As a consequence, the requirements to apply nonequilibrium work theorems are only approximately met. Starting from these considerations, we have developed a simulation scheme, called "elastic barrier dynamical freezing", according to which a stiff potential-energy barrier is enforced at the boundaries of the mobility region, preventing the particles from leaving this region of space during the nonequilibrium trajectories. The method, tested on the calculation of the distance-dependent free energy of a dimer immersed into a Lennard-Jones fluid, provides an accuracy comparable to the conventional steered molecular dynamics, with a computational speedup exceeding a few orders of magnitude.

Nonequilibrium Candidate Monte Carlo Simulations with Configurational Freezing Schemes

Nonequilibrium Candidate Monte Carlo simulation [Nilmeier et al., Proc. Natl. Acad. Sci. U.S.A. 2011, 108, E1009-E1018] is a tool devised to design Monte Carlo moves with high acceptance probabilities that connect uncorrelated configurations. Such moves are generated through nonequilibrium driven dynamics, producing candidate configurations accepted with a Monte Carlo-like criterion that preserves the equilibrium distribution. The probability of accepting a candidate configuration as the next sample in the Markov chain basically depends on the work performed on the system during the nonequilibrium trajectory and increases with decreasing such a work. It is thus strategically relevant to find ways of producing nonequilibrium moves with low work, namely moves where dissipation is as low as possible. This is the goal of our methodology, in which we combine Nonequilibrium Candidate Monte Carlo with Configurational Freezing schemes developed by Nicolini et al. (J. Chem. Theory Comput. 2011, 7, 582-593). The idea is to limit the configurational sampling to particles of a well-established region of the simulation sample, namely the region where dissipation occurs, while leaving fixed the other particles. This allows to make the system relaxation faster around the region perturbed by the finite-time switching move and hence to reduce the dissipated work, eventually enhancing the probability of accepting the generated move. Our combined approach enhances significantly configurational sampling, as shown by the case of a bistable dimer immersed in a dense fluid.

Degenerate optimal paths in thermally isolated systems

We present an analysis of the work performed on a system of interest that is kept thermally isolated during the switching of a control parameter. We show that there exists, for a certain class of systems, a finite-time family of switching protocols for which the work is equal to the quasistatic value. These optimal paths are obtained within linear response for systems initially prepared in a canonical distribution. According to our approach, such protocols are composed of a linear part plus a function which is odd with respect to time reversal. For systems with one degree of freedom, we claim that these optimal paths may also lead to the conservation of the corresponding adiabatic invariant. This points to an interesting connection between work and the conservation of the volume enclosed by the energy shell. To illustrate our findings, we solve analytically the harmonic oscillator and present numerical results for certain anharmonic examples.

Towards bulk thermodynamics via non-equilibrium methods: gaseous methane as a case study

The equation of state of bulk materials is achieved via thermodynamic derivatives of the free energy yielded by nonequilibrium transformations and Jarzynski equality.

Identifying MurI uncompetitive inhibitors by correlating decomposed binding energies with bioactivity

MurI uncompetitive inhibitors can be virtually identified by a new method that correlates decomposed binding free energies with the bioactivity.

High-precision work distributions for extreme nonequilibrium processes in large systems

The distributions of work for strongly nonequilibrium processes are studied using a very general form of a large-deviation approach, which allows one to study distributions down to extremely small probabilities of almost arbitrary quantities of interest for equilibrium, nonequilibrium stationary, and even nonstationary processes. The method is applied to quickly vary the external field in a wide range B = 3 ↔ 0 for a critical (T = 2.269) two-dimensional Ising system of size L × L = 128 × 128. To obtain free-energy differences from the work distributions, they must be studied in ranges where the probabilities are as small as 10^{-240}, which is not possible using direct simulation approaches. By comparison with the exact free energies, which are available for this model for the zero-field case, one sees that the present approach allows one to obtain the free energy with a very high relative precision of 10^{-4}. This works well also for a nonzero field, i.e., for a case where standard umbrella-sampling methods are not efficient to calculate free energies. Furthermore, for the present case it is verified that the resulting distributions of work for forward and backward processes fulfill Crooks theorem with high precision. Finally, the free energy for the Ising magnet as a function of the field strength is obtained.

Path-breaking schemes for nonequilibrium free energy calculations

We propose a path-breaking route to the enhancement of unidirectional nonequilibrium simulations for the calculation of free energy differences via Jarzynski's equality [C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997)]. One of the most important limitations of unidirectional nonequilibrium simulations is the amount of realizations necessary to reach suitable convergence of the work exponential average featuring the Jarzynski's relationship. In this respect, a significant improvement of the performances could be obtained by finding a way of stopping trajectories with negligible contribution to the work exponential average, before their normal end. This is achieved using path-breaking schemes which are essentially based on periodic checks of the work dissipated during the pulling trajectories. Such schemes can be based either on breaking trajectories whose dissipated work exceeds a given threshold or on breaking trajectories with a probability increasing with the dissipated work. In both cases, the computer time needed to carry out a series of nonequilibrium trajectories is reduced up to a factor ranging from 2 to more than 10, at least for the processes under consideration in the present study. The efficiency depends on several aspects, such as the type of process, the number of check-points along the pathway and the pulling rate as well. The method is illustrated through radically different processes, i.e., the helix-coil transition of deca-alanine and the pulling of the distance between two methane molecules in water solution.

Toward quantitative estimates of binding affinities for protein-ligand systems involving large inhibitor compounds: a steered molecular dynamics simulation route

Understanding binding mechanisms between enzymes and potential inhibitors and quantifying protein-ligand affinities in terms of binding free energy is of primary importance in drug design studies. In this respect, several approaches based on molecular dynamics simulations, often combined with docking techniques, have been exploited to investigate the physicochemical properties of complexes of pharmaceutical interest. Even if the geometric properties of a modeled protein-ligand complex can be well predicted by computational methods, it is still challenging to rank with chemical accuracy a series of ligand analogues in a consistent way. In this article, we face this issue calculating relative binding free energies of a focal adhesion kinase, an important target for the development of anticancer drugs, with pyrrolopyrimidine-based ligands having different inhibitory power. To this aim, we employ steered molecular dynamics simulations combined with nonequilibrium work theorems for free energy calculations. This technique proves very powerful when a series of ligand analogues is considered, allowing one to tackle estimation of protein-ligand relative binding free energies in a reasonable time. In our cases, the calculated binding affinities are comparable with those recovered from experiments by exploiting the Michaelis-Menten mechanism with a competitive inhibitor.

Pulling-spring modulation as a method for improving the potential-of-mean-force reconstruction in single-molecule manipulation experiments

The free-energy calculation is usually limited to close to equilibrium perturbation regimes because faster perturbations introduce a bias into the estimate. The Jarzynski equality offers a solution to this problem by directly connecting the free-energy difference and the external work, regardless how far from equilibrium that work may be. However, a limited sampling coupled to the fast perturbation introduces a slowly converging bias into the Jarzynski free-energy estimate also. In this paper we present two perturbation protocols devised with the intention to overcome the convergence issues of the Jarzynski-based potential of mean force estimation in the single-molecule, constant velocity manipulation experiments. The protocols are designed to improve the convergence issues by increasing the variation of the external work through the modulation of the spring used to pull a molecule. Of the two methods, the one which continuously changes the amplitude of the spring stiffness offers an excellent reconstruction and requires less than one tenth of the samples required by the normal, constant spring pulling to produce the same quality of the reconstruction.

Stochastic thermodynamics, fluctuation theorems and molecular machines

Stochastic thermodynamics as reviewed here systematically provides a framework for extending the notions of classical thermodynamics such as work, heat and entropy production to the level of individual trajectories of well-defined non-equilibrium ensembles. It applies whenever a non-equilibrium process is still coupled to one (or several) heat bath(s) of constant temperature. Paradigmatic systems are single colloidal particles in time-dependent laser traps, polymers in external flow, enzymes and molecular motors in single molecule assays, small biochemical networks and thermoelectric devices involving single electron transport. For such systems, a first-law like energy balance can be identified along fluctuating trajectories. For a basic Markovian dynamics implemented either on the continuum level with Langevin equations or on a discrete set of states as a master equation, thermodynamic consistency imposes a local-detailed balance constraint on noise and rates, respectively. Various integral and detailed fluctuation theorems, which are derived here in a unifying approach from one master theorem, constrain the probability distributions for work, heat and entropy production depending on the nature of the system and the choice of non-equilibrium conditions. For non-equilibrium steady states, particularly strong results hold like a generalized fluctuation-dissipation theorem involving entropy production. Ramifications and applications of these concepts include optimal driving between specified states in finite time, the role of measurement-based feedback processes and the relation between dissipation and irreversibility. Efficiency and, in particular, efficiency at maximum power can be discussed systematically beyond the linear response regime for two classes of molecular machines, isothermal ones such as molecular motors, and heat engines such as thermoelectric devices, using a common framework based on a cycle decomposition of entropy production.

Adaptive steered molecular dynamics: validation of the selection criterion and benchmarking energetics in vacuum

The potential of mean force (PMF) for stretching decaalanine in vacuum was determined earlier by Park and Schulten [J. Chem. Phys. 120, 5946 (2004)] in a landmark article demonstrating the efficacy of combining steered molecular dynamics and Jarzynski's nonequilibrium relation. In this study, the recently developed adaptive steered molecular dynamics (ASMD) algorithm [G. Ozer, E. Valeev, S. Quirk, and R. Hernandez, J. Chem. Theory Comput. 6, 3026 (2010)] is used to reproduce the PMF of the unraveling of decaalanine in vacuum by averaging over fewer nonequilibrium trajectories. The efficiency and accuracy of the method are demonstrated through the agreement with the earlier work by Park and Schulten, a series of convergence checks compared to alternate SMD pulling strategies, and an analytical proof. The nonequilibrium trajectories obtained through ASMD have also been used to analyze the intrapeptide hydrogen bonds along the stretching coordinate. As the decaalanine helix is stretched, the initially stabilized i → i + 4 contacts (α-helix) is replaced by i → i + 3 contacts (3(10)-helix). No significant formation of i → i + 5 hydrogen bonds (π-helix) is observed.

Hysteresis and nonequilibrium work theorem for DNA unzipping

We study by using Monte Carlo simulations the hysteresis in unzipping and rezipping of a double stranded DNA (dsDNA) by pulling its strands in opposite directions in the fixed force ensemble. The force is increased at a constant rate from an initial value g(0) to some maximum value g(m) that lies above the phase boundary and then decreased back again to g(0). We observed hysteresis during a complete cycle of unzipping and rezipping. We obtained probability distributions of work performed over a cycle of unzipping and rezipping for various pulling rates. The mean of the distribution is found to be close (the difference being within 10%, except for very fast pulling) to the area of the hysteresis loop. We extract the equilibrium force versus separation isotherm by using the work theorem on repeated nonequilibrium force measurements. Our method is capable of reproducing the equilibrium and the nonequilibrium force-separation isotherms for the spontaneous rezipping of dsDNA.

A Benchmark Test Set for Alchemical Free Energy Transformations and Its Use to Quantify Error in Common Free Energy Methods

There is a significant need for improved tools to validate thermophysical quantities computed via molecular simulation. In this paper we present the initial version of a benchmark set of testing methods for calculating free energies of molecular transformation in solution. This set is based on molecular changes common to many molecular design problems, such as insertion and deletion of atomic sites and changing atomic partial charges. We use this benchmark set to compare the statistical efficiency, reliability, and quality of uncertainty estimates for a number of published free energy methods, including thermodynamic integration, free energy perturbation, the Bennett acceptance ratio (BAR) and its multistate equivalent MBAR. We identify MBAR as the consistently best performing method, though other methods are frequently comparable in reliability and accuracy in many cases. We demonstrate that assumptions of Gaussian distributed errors in free energies are usually valid for most methods studied. We demonstrate that bootstrap error estimation is a robust and useful technique for estimating statistical variance for all free energy methods studied. This benchmark set is provided in a number of different file formats with the hope of becoming a useful and general tool for method comparisons.

Efficient free-energy-profile reconstruction using adaptive stochastic perturbation protocols

The Jarzynski-relation-based free-energy calculation is limited by the very slow convergence of the estimate when dissipation is high. We present two novel perturbation protocols able to significantly improve the quality of the potential of mean force (PMF) calculation by reducing the estimate's bias without increasing the number of samples. The protocols are directly applicable with both numerical simulations and real-life experiments. The improvement is achieved through the intentional but controlled widening of the distribution of the external work used to perturb a given system. Our protocols can achieve the same accuracy in PMF estimation as the normal constant velocity pulling with less than 10% of the required samples.

An Information Theory Approach to Nonlinear, Nonequilibrium Thermodynamics

Using the problem of ion channel thermodynamics as an example, we illustrate the idea of building up complex thermodynamic models by successively adding physical information. We present a new formulation of information algebra that generalizes methods of both information theory and statistical mechanics. From this foundation we derive a theory for ion channel kinetics, identifying a nonequilibrium 'process' free energy functional in addition to the well-known integrated work functionals. The Gibbs-Maxwell relation for the free energy functional is a Green-Kubo relation, applicable arbitrarily far from equilibrium, that captures the effect of non-local and time-dependent behavior from transient thermal and mechanical driving forces. Comparing the physical significance of the Lagrange multipliers to the canonical ensemble suggests definitions of nonequilibrium ensembles at constant capacitance or inductance in addition to constant resistance. Our result is that statistical mechanical descriptions derived from a few primitive algebraic operations on information can be used to create experimentally-relevant and computable models. By construction, these models may use information from more detailed atomistic simulations. Two surprising consequences to be explored in further work are that (in)distinguishability factors are automatically predicted from the problem formulation and that a direct analogue of the second law for thermodynamic entropy production is found by considering information loss in stochastic processes. The information loss identifies a novel contribution from the instantaneous information entropy that ensures non-negative loss.

Use of Umbrella Sampling to Calculate the Entrance/Exit Pathway for Z-Pro-Prolinal Inhibitor in Prolyl Oligopeptidase

Prolyl oligopeptidase (POP), a member of the prolyl endopeptidase family, is known to play a role in several neurological disorders. Its primary function is to cleave a wide range of small oligopeptides, including neuroactive peptides. We have used force biased molecular dynamics simulation to study the binding mechanism of POP. We examined three possible binding pathways using Steered Molecular Dynamics (SMD) and Umbrella Sampling (US) on a crystal structure of porcine POP with bound Z-pro-prolinal (ZPP). Using SMD, an exit pathway between the first and seventh blade of the β-propeller domain of POP was found to be a nonviable route. US on binding pathways through the β-propeller tunnel and the TYR190-GLN208 flexible loop at the interface between both POP domains allowed us to isolate the flexible loop pathway as the most probable. Further analysis of that pathway suggests a long-range covariation of the interdomain H-bond network, which indicates the possibility of large-scale domain reorientation observed in bacterial homologues and hypothesized to also occur in human POP.

Measuring the convergence of Monte Carlo free-energy calculations

The nonequilibrium work fluctuation theorem provides the way for calculations of (equilibrium) free-energy based on work measurements of nonequilibrium, finite-time processes, and their reversed counterparts by applying Bennett's acceptance ratio method. A nice property of this method is that each free-energy estimate readily yields an estimate of the asymptotic mean square error. Assuming convergence, it is easy to specify the uncertainty of the results. However, sample sizes have often to be balanced with respect to experimental or computational limitations and the question arises whether available samples of work values are sufficiently large in order to ensure convergence. Here, we propose a convergence measure for the two-sided free-energy estimator and characterize some of its properties, explain how it works, and test its statistical behavior. In total, we derive a convergence criterion for Bennett's acceptance ratio method.

Optimum protocol for fast-switching free-energy calculations

Free-energy differences computed from fast-switching simulations or measurements according to the Jarzynski equation are independent of the particular protocol specifying how the control parameter is changed in time. In contrast, the average work carried out on the system as well the accuracy of the resulting free energy strongly depend on the protocol. Recently, Schmiedl and Seifert [Phys. Rev. Lett. 98, 108301 (2007)] found that protocols that minimize the average work for a given duration of the switching process have discrete steps at the beginning and the end. Here we determine numerically the protocols that minimize the statistical error in the free energy estimate and find that such minimum error protocols have similar discrete jumps. Our analysis shows that the reduction in computational effort achieved by the use of steplike protocols can be considerable. Such large savings of computing time, however, typically occur for parameter ranges in which an application of the Jarzynski equation is impractical due to large statistical errors arising from the exponential work average.

Optimal estimators and asymptotic variances for nonequilibrium path-ensemble averages

Existing optimal estimators of nonequilibrium path-ensemble averages are shown to fall within the framework of extended bridge sampling. Using this framework, we derive a general minimal-variance estimator that can combine nonequilibrium trajectory data sampled from multiple path-ensembles to estimate arbitrary functions of nonequilibrium expectations. The framework is also applied to obtain asymptotic variance estimates, which are a useful measure of statistical uncertainty. In particular, we develop asymptotic variance estimates pertaining to Jarzynski's equality for free energies and the Hummer-Szabo expressions for the potential of mean force, calculated from uni- or bidirectional path samples. These estimators are demonstrated on a model single-molecule pulling experiment. In these simulations, the asymptotic variance expression is found to accurately characterize the confidence intervals around estimators when the bias is small. Hence, the confidence intervals are inaccurately described for unidirectional estimates with large bias, but for this model it largely reflects the true error in a bidirectional estimator derived by Minh and Adib.

Improving fast-switching free energy estimates by dynamical freezing

An important limitation of nonequilibrium pulling experiments/simulations in recovering free energy differences is the poor convergence of path-ensemble averages. Therefore, a large number of fast-switching trajectories needs to achieve free energy estimates with acceptable accuracy. We propose a method to improve free energy estimates by drastically lowering the computational cost of steered molecular dynamics simulations employed to realize such trajectories. This is accomplished by generating trajectories where the particles not directly involved in the driven process are dynamically frozen. Such a freezing is dynamical rather than thermal because it is reached by a synchronous scaling of atomic masses and velocities keeping the kinetic energy of each particle unchanged. The forces between dynamically frozen particles can then be calculated rarely. Thus, it is possible to generate realizations of a process whose computational cost is not correlated with the size of the whole system, but only with that of the reaction site. The method is illustrated on a simple model system and its general applicability is discussed.

Characteristic of Bennett's acceptance ratio method

A powerful and well-established tool for free-energy estimation is Bennett's acceptance ratio method. Central properties of this estimator, which employs samples of work values of a forward and its time-reversed process, are known: for given sets of measured work values, it results in the best estimate of the free-energy difference in the large sample limit. Here we state and prove a further characteristic of the acceptance ratio method: the convexity of its mean-square error. As a two-sided estimator, it depends on the ratio of the numbers of forward and reverse work values used. Convexity of its mean-square error immediately implies that there exists a unique optimal ratio for which the error becomes minimal. Further, it yields insight into the relation of the acceptance ratio method and estimators based on the Jarzynski equation. As an application, we study the performance of a dynamic strategy of sampling forward and reverse work values.

Asymptotics of work distributions in nonequilibrium systems

The asymptotic behavior of the work distribution in driven nonequilibrium systems is determined using the method of optimal fluctuations. For systems described by Langevin dynamics the corresponding Euler-Lagrange equation together with the appropriate boundary conditions and an equation for the leading pre-exponential factor are derived. The method is applied to three representative examples and the results are used to improve the accuracy of free-energy estimates based on the application of the Jarzynski equation.

Density-dependent analysis of nonequilibrium paths improves free energy estimates

When a system is driven out of equilibrium by a time-dependent protocol that modifies the Hamiltonian, it follows a nonequilibrium path. Samples of these paths can be used in nonequilibrium work theorems to estimate equilibrium quantities such as free energy differences. Here, we consider analyzing paths generated with one protocol using another one. It is posited that analysis protocols which minimize the lag, the difference between the nonequilibrium and the instantaneous equilibrium densities, will reduce the dissipation of reprocessed trajectories and lead to better free energy estimates. Indeed, when minimal lag analysis protocols based on exactly soluble propagators or relative entropies are applied to several test cases, substantial gains in the accuracy and precision of estimated free energy differences are observed.

Velocity convergence of free energy surfaces from single-molecule measurements using Jarzynski's equality

We studied the velocity dependence of mechanical unfolding of single protein molecules with the atomic force microscope. We showed that with enough realizations, the free energy surfaces reconstructed from Jarzynski's equality converge with respect to pulling velocity, in good agreement with theory. Using the I27 domain of titin as an example, we estimated the required number of realizations for a given pulling velocity, and suggested the optimal range of velocities for single-molecule experiments. The results demonstrate that Jarzynski's equality is a powerful and practical tool for reconstructing free energy landscapes.

Optimal protocols for minimal work processes in underdamped stochastic thermodynamics

For systems in an externally controllable time-dependent potential, the optimal protocol minimizes the mean work spent in a finite-time transition between two given equilibrium states. For overdamped dynamics which ignores inertia effects, the optimal protocol has been found to involve jumps of the control parameter at the beginning and end of the process. Including the inertia term, we show that this feature not only persists but that even delta-peak-like changes of the control parameter at both boundaries make the process optimal. These results are obtained by analyzing two simple paradigmatic cases: First, a Brownian particle dragged by a harmonic optical trap through a viscous fluid and, second, a Brownian particle subject to an optical trap with time-dependent stiffness. These insights could be used to improve free energy calculations via either thermodynamic integration or "fast growth" methods using Jarzynski's equality.

Protein-protein interaction investigated by steered molecular dynamics: the TCR-pMHC complex

We present a novel steered molecular dynamics scheme to induce the dissociation of large protein-protein complexes. We apply this scheme to study the interaction of a T cell receptor (TCR) with a major histocompatibility complex (MHC) presenting a peptide (p). Two TCR-pMHC complexes are considered, which only differ by the mutation of a single amino acid on the peptide; one is a strong agonist that produces T cell activation in vivo, while the other is an antagonist. We investigate the interaction mechanism from a large number of unbinding trajectories by analyzing van der Waals and electrostatic interactions and by computing energy changes in proteins and solvent. In addition, dissociation potentials of mean force are calculated with the Jarzynski identity, using an averaging method developed for our steering scheme. We analyze the convergence of the Jarzynski exponential average, which is hampered by the large amount of dissipative work involved and the complexity of the system. The resulting dissociation free energies largely underestimate experimental values, but the simulations are able to clearly differentiate between wild-type and mutated TCR-pMHC and give insights into the dissociation mechanism.

Escorted free energy simulations: improving convergence by reducing dissipation

Nonequilibrium, "fast switching" estimates of equilibrium free energy differences DeltaF are often plagued by poor convergence due to dissipation. We propose a method to improve these estimates by generating trajectories with reduced dissipation. Introducing an artificial flow field that couples the system coordinates to the external parameter driving the simulation, we derive an identity for DeltaF in terms of the resulting trajectories. When the flow field effectively escorts the system along a near-equilibrium path, the free energy estimate converges efficiently and accurately. We illustrate our method on a model system and discuss the general applicability of our approach.

Optimized free energies from bidirectional single-molecule force spectroscopy

An optimized method for estimating path-ensemble averages using data from processes driven in opposite directions is presented. Based on this estimator, bidirectional expressions for reconstructing free energies and potentials of mean force from single-molecule force spectroscopy-valid for biasing potentials of arbitrary stiffness-are developed. Numerical simulations on a model potential indicate that these methods perform better than unidirectional strategies.