Protein-Ligand Binding Free-Energy Calculations with ARROW─A Purely First-Principles Parameterized Polarizable Force Field.

Protein-ligand binding free-energy calculations using molecular dynamics (MD) simulations have emerged as a powerful tool for in silico drug design. Here, we present results obtained with the ARROW force field (FF)─a multipolar polarizable and physics-based model with all parameters fitted entirely to high-level ab initio quantum mechanical (QM) calculations. ARROW has already proven its ability to determine solvation free energy of arbitrary neutral compounds with unprecedented accuracy. The ARROW FF parameterization is now extended to include coverage of all amino acids including charged groups, allowing molecular simulations of a series of protein-ligand systems and prediction of their relative binding free energies. We ensure adequate sampling by applying a novel technique that is based on coupling the Hamiltonian Replica exchange (HREX) with a conformation reservoir generated via potential softening and nonequilibrium MD. ARROW provides predictions with near chemical accuracy (mean absolute error of ∼0.5 kcal/mol) for two of the three protein systems studied here (MCL1 and Thrombin). The third protein system (CDK2) reveals the difficulty in accurately describing dimer interaction energies involving polar and charged species. Overall, for all of the three protein systems studied here, ARROW FF predicts relative binding free energies of ligands with a similar accuracy level as leading nonpolarizable force fields.

Quantifying Configuration-Sampling Error in Langevin Simulations of Complex Molecular Systems.

While Langevin integrators are popular in the study of equilibrium properties of complex systems, it is challenging to estimate the timestep-induced discretization error: the degree to which the sampled phase-space or configuration-space probability density departs from the desired target density due to the use of a finite integration timestep. Sivak et al., introduced a convenient approach to approximating a natural measure of error between the sampled density and the target equilibrium density, the Kullback-Leibler (KL) divergence, in phase space, but did not specifically address the issue of configuration-space properties, which are much more commonly of interest in molecular simulations. Here, we introduce a variant of this near-equilibrium estimator capable of measuring the error in the configuration-space marginal density, validating it against a complex but exact nested Monte Carlo estimator to show that it reproduces the KL divergence with high fidelity. To illustrate its utility, we employ this new near-equilibrium estimator to assess a claim that a recently proposed Langevin integrator introduces extremely small configuration-space density errors up to the stability limit at no extra computational expense. Finally, we show how this approach to quantifying sampling bias can be applied to a wide variety of stochastic integrators by following a straightforward procedure to compute the appropriate shadow work, and describe how it can be extended to quantify the error in arbitrary marginal or conditional distributions of interest.

Geometric approach to optimal nonequilibrium control: Minimizing dissipation in nanomagnetic spin systems.

Optimal control of nanomagnets has become an urgent problem for the field of spintronics as technological tools approach thermodynamically determined limits of efficiency. In complex, fluctuating systems, such as nanomagnetic bits, finding optimal protocols is challenging, requiring detailed information about the dynamical fluctuations of the controlled system. We provide a physically transparent derivation of a metric tensor for which the length of a protocol is proportional to its dissipation. This perspective simplifies nonequilibrium optimization problems by recasting them in a geometric language. We then describe a numerical method, an instance of geometric minimum action methods, that enables computation of geodesics even when the number of control parameters is large. We apply these methods to two models of nanomagnetic bits: a Landau-Lifshitz-Gilbert description of a single magnetic spin controlled by two orthogonal magnetic fields, and a two-dimensional Ising model in which the field is spatially controlled. These calculations reveal nontrivial protocols for bit erasure and reversal, providing important, experimentally testable predictions for ultra-low-power computing.

Thermodynamic geometry of minimum-dissipation driven barrier crossing.

We explore the thermodynamic geometry of a simple system that models the bistable dynamics of nucleic acid hairpins in single molecule force-extension experiments. Near equilibrium, optimal (minimum-dissipation) driving protocols are governed by a generalized linear response friction coefficient. Our analysis demonstrates that the friction coefficient of the driving protocols is sharply peaked at the interface between metastable regions, which leads to minimum-dissipation protocols that drive rapidly within a metastable basin, but then linger longest at the interface, giving thermal fluctuations maximal time to kick the system over the barrier. Intuitively, the same principle applies generically in free energy estimation (both in steered molecular dynamics simulations and in single-molecule experiments), provides a design principle for the construction of thermodynamically efficient coupling between stochastic objects, and makes a prediction regarding the construction of evolved biomolecular motors.

Near-optimal protocols in complex nonequilibrium transformations.

The development of sophisticated experimental means to control nanoscale systems has motivated efforts to design driving protocols that minimize the energy dissipated to the environment. Computational models are a crucial tool in this practical challenge. We describe a general method for sampling an ensemble of finite-time, nonequilibrium protocols biased toward a low average dissipation. We show that this scheme can be carried out very efficiently in several limiting cases. As an application, we sample the ensemble of low-dissipation protocols that invert the magnetization of a 2D Ising model and explore how the diversity of the protocols varies in response to constraints on the average dissipation. In this example, we find that there is a large set of protocols with average dissipation close to the optimal value, which we argue is a general phenomenon.

Optimal control in nonequilibrium systems: Dynamic Riemannian geometry of the Ising model.

A general understanding of optimal control in nonequilibrium systems would illuminate the operational principles of biological and artificial nanoscale machines. Recent work has shown that a system driven out of equilibrium by a linear response protocol is endowed with a Riemannian metric related to generalized susceptibilities, and that geodesics on this manifold are the nonequilibrium control protocols with the lowest achievable dissipation. While this elegant mathematical framework has inspired numerous studies of exactly solvable systems, no description of the thermodynamic geometry yet exists when the metric cannot be derived analytically. Herein, we numerically construct the dynamic metric of the two-dimensional Ising model in order to study optimal protocols for reversing the net magnetization.

Scaling laws governing stochastic growth and division of single bacterial cells.

Uncovering the quantitative laws that govern the growth and division of single cells remains a major challenge. Using a unique combination of technologies that yields unprecedented statistical precision, we find that the sizes of individual Caulobacter crescentus cells increase exponentially in time. We also establish that they divide upon reaching a critical multiple (≈ 1.8) of their initial sizes, rather than an absolute size. We show that when the temperature is varied, the growth and division timescales scale proportionally with each other over the physiological temperature range. Strikingly, the cell-size and division-time distributions can both be rescaled by their mean values such that the condition-specific distributions collapse to universal curves. We account for these observations with a minimal stochastic model that is based on an autocatalytic cycle. It predicts the scalings, as well as specific functional forms for the universal curves. Our experimental and theoretical analysis reveals a simple physical principle governing these complex biological processes: a single temperature-dependent scale of cellular time governs the stochastic dynamics of growth and division in balanced growth conditions.

Universality in stochastic exponential growth.

Recent imaging data for single bacterial cells reveal that their mean sizes grow exponentially in time and that their size distributions collapse to a single curve when rescaled by their means. An analogous result holds for the division-time distributions. A model is needed to delineate the minimal requirements for these scaling behaviors. We formulate a microscopic theory of stochastic exponential growth as a Master Equation that accounts for these observations, in contrast to existing quantitative models of stochastic exponential growth (e.g., the Black-Scholes equation or geometric Brownian motion). Our model, the stochastic Hinshelwood cycle (SHC), is an autocatalytic reaction cycle in which each molecular species catalyzes the production of the next. By finding exact analytical solutions to the SHC and the corresponding first passage time problem, we uncover universal signatures of fluctuations in exponential growth and division. The model makes minimal assumptions, and we describe how more complex reaction networks can reduce to such a cycle. We thus expect similar scalings to be discovered in stochastic processes resulting in exponential growth that appear in diverse contexts such as cosmology, finance, technology, and population growth.

Time step rescaling recovers continuous-time dynamical properties for discrete-time Langevin integration of nonequilibrium systems.

When simulating molecular systems using deterministic equations of motion (e.g., Newtonian dynamics), such equations are generally numerically integrated according to a well-developed set of algorithms that share commonly agreed-upon desirable properties. However, for stochastic equations of motion (e.g., Langevin dynamics), there is still broad disagreement over which integration algorithms are most appropriate. While multiple desiderata have been proposed throughout the literature, consensus on which criteria are important is absent, and no published integration scheme satisfies all desiderata simultaneously. Additional nontrivial complications stem from simulating systems driven out of equilibrium using existing stochastic integration schemes in conjunction with recently developed nonequilibrium fluctuation theorems. Here, we examine a family of discrete time integration schemes for Langevin dynamics, assessing how each member satisfies a variety of desiderata that have been enumerated in prior efforts to construct suitable Langevin integrators. We show that the incorporation of a novel time step rescaling in the deterministic updates of position and velocity can correct a number of dynamical defects in these integrators. Finally, we identify a particular splitting (related to the velocity Verlet discretization) that has essentially universally appropriate properties for the simulation of Langevin dynamics for molecular systems in equilibrium, nonequilibrium, and path sampling contexts.

Geometry of thermodynamic control.

A deeper understanding of nonequilibrium phenomena is needed to reveal the principles governing natural and synthetic molecular machines. Recent work has shown that when a thermodynamic system is driven from equilibrium then, in the linear response regime, the space of controllable parameters has a Riemannian geometry induced by a generalized friction tensor. We exploit this geometric insight to construct closed-form expressions for minimal-dissipation protocols for a particle diffusing in a one-dimensional harmonic potential, where the spring constant, inverse temperature, and trap location are adjusted simultaneously. These optimal protocols are geodesics on the Riemannian manifold and reveal that this simple model has a surprisingly rich geometry. We test these optimal protocols via a numerical implementation of the Fokker-Planck equation and demonstrate that the friction tensor arises naturally from a first-order expansion in temporal derivatives of the control parameters, without appealing directly to linear response theory.

Thermodynamic metrics and optimal paths.

A fundamental problem in modern thermodynamics is how a molecular-scale machine performs useful work, while operating away from thermal equilibrium without excessive dissipation. To this end, we derive a friction tensor that induces a Riemannian manifold on the space of thermodynamic states. Within the linear-response regime, this metric structure controls the dissipation of finite-time transformations, and bestows optimal protocols with many useful properties. We discuss the connection to the existing thermodynamic length formalism, and demonstrate the utility of this metric by solving for optimal control parameter protocols in a simple nonequilibrium model.

Thermodynamics of prediction.

A system responding to a stochastic driving signal can be interpreted as computing, by means of its dynamics, an implicit model of the environmental variables. The system's state retains information about past environmental fluctuations, and a fraction of this information is predictive of future ones. The remaining nonpredictive information reflects model complexity that does not improve predictive power, and thus represents the ineffectiveness of the model. We expose the fundamental equivalence between this model inefficiency and thermodynamic inefficiency, measured by dissipation. Our results hold arbitrarily far from thermodynamic equilibrium and are applicable to a wide range of systems, including biomolecular machines. They highlight a profound connection between the effective use of information and efficient thermodynamic operation: any system constructed to keep memory about its environment and to operate with maximal energetic efficiency has to be predictive.

Near-equilibrium measurements of nonequilibrium free energy.

A central endeavor of thermodynamics is the measurement of free energy changes. Regrettably, although we can measure the free energy of a system in thermodynamic equilibrium, typically all we can say about the free energy of a nonequilibrium ensemble is that it is larger than that of the same system at equilibrium. Herein, we derive a formally exact expression for the probability distribution of a driven system, which involves path ensemble averages of the work over trajectories of the time-reversed system. From this we find a simple near-equilibrium approximation for the free energy in terms of an excess mean time-reversed work, which can be experimentally measured on real systems. With analysis and computer simulation, we demonstrate the accuracy of our approximations for several simple models.

Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation.

Metropolis Monte Carlo simulation is a powerful tool for studying the equilibrium properties of matter. In complex condensed-phase systems, however, it is difficult to design Monte Carlo moves with high acceptance probabilities that also rapidly sample uncorrelated configurations. Here, we introduce a new class of moves based on nonequilibrium dynamics: Candidate configurations are generated through a finite-time process in which a system is actively driven out of equilibrium, and accepted with criteria that preserve the equilibrium distribution. The acceptance rule is similar to the Metropolis acceptance probability, but related to the nonequilibrium work rather than the instantaneous energy difference. Our method is applicable to sampling from both a single thermodynamic state or a mixture of thermodynamic states, and allows both coordinates and thermodynamic parameters to be driven in nonequilibrium proposals. Whereas generating finite-time switching trajectories incurs an additional cost, driving some degrees of freedom while allowing others to evolve naturally can lead to large enhancements in acceptance probabilities, greatly reducing structural correlation times. Using nonequilibrium driven processes vastly expands the repertoire of useful Monte Carlo proposals in simulations of dense solvated systems.

Self-organization of the Escherichia coli chemotaxis network imaged with super-resolution light microscopy.

The Escherichia coli chemotaxis network is a model system for biological signal processing. In E. coli, transmembrane receptors responsible for signal transduction assemble into large clusters containing several thousand proteins. These sensory clusters have been observed at cell poles and future division sites. Despite extensive study, it remains unclear how chemotaxis clusters form, what controls cluster size and density, and how the cellular location of clusters is robustly maintained in growing and dividing cells. Here, we use photoactivated localization microscopy (PALM) to map the cellular locations of three proteins central to bacterial chemotaxis (the Tar receptor, CheY, and CheW) with a precision of 15 nm. We find that cluster sizes are approximately exponentially distributed, with no characteristic cluster size. One-third of Tar receptors are part of smaller lateral clusters and not of the large polar clusters. Analysis of the relative cellular locations of 1.1 million individual proteins (from 326 cells) suggests that clusters form via stochastic self-assembly. The super-resolution PALM maps of E. coli receptors support the notion that stochastic self-assembly can create and maintain approximately periodic structures in biological membranes, without direct cytoskeletal involvement or active transport.

Far-from-equilibrium measurements of thermodynamic length.

Thermodynamic length is a path function that generalizes the notion of length to the surface of thermodynamic states. Here, we show how to measure thermodynamic length in far-from-equilibrium experiments using the work fluctuation relations. For these microscopic systems, it proves necessary to define the thermodynamic length in terms of the Fisher information. Consequently, the thermodynamic length can be directly related to the magnitude of fluctuations about equilibrium. The work fluctuation relations link the work and the free-energy change during an external perturbation on a system. We use this result to determine equilibrium averages at intermediate points of the protocol in which the system is out of equilibrium. This allows us to extend Bennett's method to determine the potential of the mean force, as well as the thermodynamic length, in single-molecule experiments.

Length of time's arrow.

An unresolved problem in physics is how the thermodynamic arrow of time arises from an underlying time reversible dynamics. We contribute to this issue by developing a measure of time-symmetry breaking, and by using the work fluctuation relations, we determine the time asymmetry of recent single molecule RNA unfolding experiments. We define time asymmetry as the Jensen-Shannon divergence between trajectory probability distributions of an experiment and its time-reversed conjugate. Among other interesting properties, the length of time's arrow bounds the average dissipation and determines the difficulty of accurately estimating free energy differences in nonequilibrium experiments.

Bayesian estimates of free energies from nonequilibrium work data in the presence of instrument noise.

The Jarzynski equality and the fluctuation theorem relate equilibrium free energy differences to nonequilibrium measurements of the work. These relations extend to single-molecule experiments that have probed the finite-time thermodynamics of proteins and nucleic acids. The effects of experimental error and instrument noise have not been considered previously. Here, we present a Bayesian formalism for estimating free energy changes from nonequilibrium work measurements that compensates for instrument noise and combines data from multiple driving protocols. We reanalyze a recent set of experiments in which a single RNA hairpin is unfolded and refolded using optical tweezers at three different rates. Interestingly, the fastest and farthest-from-equilibrium measurements contain the least instrumental noise and, therefore, provide a more accurate estimate of the free energies than a few slow, more noisy, near-equilibrium measurements. The methods we propose here will extend the scope of single-molecule experiments; they can be used in the analysis of data from measurements with atomic force microscopy, optical, and magnetic tweezers.

Measuring thermodynamic length.

Thermodynamic length is a metric distance between equilibrium thermodynamic states. Among other interesting properties, this metric asymptotically bounds the dissipation induced by a finite time transformation of a thermodynamic system. It is also connected to the Jensen-Shannon divergence, Fisher information, and Rao's entropy differential metric. Therefore, thermodynamic length is of central interest in understanding matter out of equilibrium. In this Letter, we will consider how to define thermodynamic length for a small system described by equilibrium statistical mechanics and how to measure thermodynamic length within a computer simulation. Surprisingly, Bennett's classic acceptance ratio method for measuring free energy differences also measures thermodynamic length.

Beyond Boltzmann-Gibbs statistics: maximum entropy hyperensembles out of equilibrium.

What is the best description that we can construct of a thermodynamic system that is not in equilibrium, given only one, or a few, extra parameters over and above those needed for a description of the same system at equilibrium? Here, we argue the most appropriate additional parameter is the nonequilibrium entropy of the system. Moreover, we should not attempt to estimate the probability distribution of the system directly, but rather the metaprobability (or hyperensemble) that the system is described by a particular probability distribution. The result is an entropic distribution with two parameters, one a nonequilibrium temperature, and the other a measure of distance from equilibrium. This dispersion parameter smoothly interpolates between certainty of a canonical distribution at equilibrium and great uncertainty as to the probability distribution as we move away from equilibrium. We deduce that, in general, large, rare fluctuations become far more common as we move away from equilibrium.