Markov chain Monte Carlo method without detailed balance

Control of probability flow in Markov chain Monte Carlo-Nonreversibility and lifting

The Markov chain Monte Carlo (MCMC) method is widely used in various fields as a powerful numerical integration technique for systems with many degrees of freedom. In MCMC methods, probabilistic state transitions can be considered as a random walk in state space, and random walks allow for sampling from complex distributions. However, paradoxically, it is necessary to carefully suppress the randomness of the random walk to improve computational efficiency. By breaking detailed balance, we can create a probability flow in the state space and perform more efficient sampling along this flow. Motivated by this idea, practical and efficient nonreversible MCMC methods have been developed over the past ten years. In particular, the lifting technique, which introduces probability flows in an extended state space, has been applied to various systems and has proven more efficient than conventional reversible updates. We review and discuss several practical approaches to implementing nonreversible MCMC methods, including the shift method in the cumulative distribution and the directed-worm algorithm.

Why Is Arginine the Only Amino Acid That Inhibits Polyglutamine Monomers from Taking on Toxic Conformations?

Polyglutamine (polyQ) diseases are devastating neurodegenerative disorders characterized by abnormal expansion of glutamine repeats within specific proteins. The aggregation of polyQ proteins is a critical pathological hallmark of these diseases. Arginine was identified as a promising inhibitory compound because it prevents polyQ-protein monomers from forming intra- and intermolecular β-sheet structures and hinders polyQ proteins from aggregating to form oligomers. Such an aggregation inhibitory effect was not observed in other amino acids. However, the underlying molecular mechanism of the aggregation inhibition and the factors that differentiate arginine from other amino acids, in terms of the inhibition of the polyQ-protein aggregation, remain poorly understood. Here, we performed replica-permutation molecular dynamics simulations to elucidate the molecular mechanism by which arginine inhibits the formation of the intramolecular β-sheet structure of a polyQ monomer. We found that the intramolecular β-sheet structure with more than four β-bridges of the polyQ monomer with arginine is more unstable than without any ligand and with lysine. We also found that arginine has 1.6-2.1 times more contact with polyQ than lysine. In addition, we revealed that arginine forms more hydrogen bonds with the main chain of the polyQ monomer than lysine. More hydrogen bonds formed between arginine and polyQ inhibit polyQ from forming the long intramolecular β-sheet structure. It is known that intramolecular β-sheet structure enhances intermolecular β-sheet structure between proteins. These effects are thought to be the reason for the inhibition of polyQ aggregation. This study provides insights into the molecular events underlying arginine's inhibition of polyQ-protein aggregation.

Advanced computational approaches to understand protein aggregation

Protein aggregation is a widespread phenomenon implicated in debilitating diseases like Alzheimer's, Parkinson's, and cataracts, presenting complex hurdles for the field of molecular biology. In this review, we explore the evolving realm of computational methods and bioinformatics tools that have revolutionized our comprehension of protein aggregation. Beginning with a discussion of the multifaceted challenges associated with understanding this process and emphasizing the critical need for precise predictive tools, we highlight how computational techniques have become indispensable for understanding protein aggregation. We focus on molecular simulations, notably molecular dynamics (MD) simulations, spanning from atomistic to coarse-grained levels, which have emerged as pivotal tools in unraveling the complex dynamics governing protein aggregation in diseases such as cataracts, Alzheimer's, and Parkinson's. MD simulations provide microscopic insights into protein interactions and the subtleties of aggregation pathways, with advanced techniques like replica exchange molecular dynamics, Metadynamics (MetaD), and umbrella sampling enhancing our understanding by probing intricate energy landscapes and transition states. We delve into specific applications of MD simulations, elucidating the chaperone mechanism underlying cataract formation using Markov state modeling and the intricate pathways and interactions driving the toxic aggregate formation in Alzheimer's and Parkinson's disease. Transitioning we highlight how computational techniques, including bioinformatics, sequence analysis, structural data, machine learning algorithms, and artificial intelligence have become indispensable for predicting protein aggregation propensity and locating aggregation-prone regions within protein sequences. Throughout our exploration, we underscore the symbiotic relationship between computational approaches and empirical data, which has paved the way for potential therapeutic strategies against protein aggregation-related diseases. In conclusion, this review offers a comprehensive overview of advanced computational methodologies and bioinformatics tools that have catalyzed breakthroughs in unraveling the molecular basis of protein aggregation, with significant implications for clinical interventions, standing at the intersection of computational biology and experimental research.

Perspective for Molecular Dynamics Simulation Studies of Amyloid-β Aggregates

The cause of Alzheimer's disease is related to aggregates such as oligomers and amyloid fibrils consisting of amyloid-β (Aβ) peptides. Molecular dynamics (MD) simulation studies have been conducted to understand the molecular mechanism of the formation and disruption of Aβ aggregates. In this Perspective, the MD simulation studies are classified into four categories, focusing on the target systems: aggregation of Aβ peptides in bulk solution, Aβ aggregation at the interface, aggregation inhibitor against Aβ peptides, and nonequilibrium MD simulation of Aβ aggregates. MD simulation studies in these categories are first reviewed. Future perspectives in each category are then presented. Finally, the overall perspective is presented on how MD simulations of Aβ aggregates can be utilized for developing Alzheimer's disease treatment.

Thermodynamic Bounds on Correlation Times

We derive a variational expression for the correlation time of physical observables in steady-state diffusive systems. As a consequence of this variational expression, we obtain lower bounds on the correlation time, which provide speed limits on the self-averaging of observables. In equilibrium, the bound takes the form of a trade-off relation between the long- and short-time fluctuations of an observable. Out of equilibrium, the trade-off can be violated, leading to an acceleration of self-averaging. We relate this violation to the steady-state entropy production rate, as well as the geometric structure of the irreversible currents, giving rise to two complementary speed limits. One of these can be formulated as a lower estimate on the entropy production from the measurement of time-symmetric observables. Using an illustrating example, we show the intricate behavior of the correlation time out of equilibrium for different classes of observables and how this can be used to partially infer dissipation even if no time-reversal symmetry breaking can be observed in the trajectories of the observable.

Ingenuity in performing replica permutation: How to order the state labels for improving sampling efficiency

In the replica-permutation method, an advanced version of the replica-exchange method, all combinations of replicas and parameters are considered for parameter permutation, and a list of all the combinations is prepared. Here, we report that the temperature transition probability depends on how the list is created, especially in replica permutation with solute tempering (RPST). We found that the transition probabilities decrease at large replica indices when the combinations are sequentially assigned to the state labels as in the originally proposed list. To solve this problem, we propose to modify the list by randomly assigning the combinations to the state labels. We performed molecular dynamics simulations of amyloid-β(16-22) peptides using RPST with the "randomly assigned" list (RPST-RA) and RPST with the "sequentially assigned" list (RPST-SA). The results show the decreases in the transition probabilities in RPST-SA are eliminated, and the sampling efficiency is improved in RPST-RA.

Lifted directed-worm algorithm

Nonreversible Markov chains can outperform reversible chains in the Markov chain Monte Carlo method. Lifting is a versatile approach to introducing net stochastic flow in state space and constructing a nonreversible Markov chain. We present here an application of the lifting technique to the directed-worm algorithm. The transition probability of the worm update is optimized using the geometric allocation approach; the worm backscattering probability is minimized, and the stochastic flow breaking the detailed balance is maximized. We demonstrate the performance improvement over the previous worm and cluster algorithms for the four-dimensional hypercubic lattice Ising model. The sampling efficiency of the present algorithm is approximately 80, 5, and 1.7 times as high as those of the standard worm algorithm, the Wolff cluster algorithm, and the previous lifted worm algorithm, respectively. We estimate the dynamic critical exponent of the hypercubic lattice Ising model to be z≈0 in the worm and the Wolff cluster updates. The lifted version of the directed-worm algorithm can be applied to a variety of quantum systems as well as classical systems.

Key Residue for Aggregation of Amyloid-β Peptides

It is known that oligomers of amyloid-β (Aβ) peptide are associated with Alzheimer's disease. Aβ has two isoforms: Aβ40 and Aβ42. Although the difference between Aβ40 and Aβ42 is only two additional C-terminal residues, Aβ42 aggregates much faster than Aβ40. It is unknown what role the C-terminal two residues play in accelerating aggregation. Since Aβ42 is more toxic than Aβ40, its oligomerization process needs to be clarified. Moreover, clarifying the differences between the oligomerization processes of Aβ40 and Aβ42 is essential to elucidate the key factors of oligomerization. Therefore, to investigate the dimerization process, which is the early oligomerization process, Hamiltonian replica-permutation molecular dynamics simulations were performed for Aβ40 and Aβ42. We identified a key residue, Arg5, for the Aβ42 dimerization. The two additional residues in Aβ42 allow the C-terminus to form contact with Arg5 because of the electrostatic attraction between them, and this contact stabilizes the β-hairpin. This β-hairpin promotes dimer formation through the intermolecular β-bridges. Thus, we examined the effects of amino acid substitutions of Arg5, thereby confirming that the mutations remarkably suppressed the aggregation of Aβ42. Moreover, the mutations of Arg5 suppressed the Aβ40 aggregation. It was found by analyzing the simulations that Arg5 is important for Aβ40 to form intermolecular contacts. Thus, it was clarified that the role of Arg5 in the oligomerization process varies due to the two additional C-terminal residues.

Accelerated Jarzynski estimator with deterministic virtual trajectories

The Jarzynski estimator is a powerful tool that uses nonequilibrium statistical physics to numerically obtain partition functions of probability distributions. The estimator reconstructs partition functions with trajectories of the simulated Langevin dynamics through the Jarzynski equality. However, the original estimator suffers from slow convergence because it depends on rare trajectories of stochastic dynamics. In this paper, we present a method to significantly accelerate the convergence by introducing deterministic virtual trajectories generated in augmented state space under the Hamiltonian dynamics. We theoretically show that our approach achieves second-order acceleration compared to a naive estimator with the Langevin dynamics and zero variance estimation on harmonic potentials. We also present numerical experiments on three multimodal distributions and a practical example in which the proposed method outperforms the conventional method, and we provide theoretical explanations.

Accelerating, to some extent, the p-spin dynamics

We consider a detailed-balance-violating dynamics whose stationary state is a prescribed Boltzmann distribution. Such dynamics have been shown to be faster than any equilibrium counterpart. We quantify the gain in convergence speed for a system whose energy landscape displays one and then an infinite number of energy barriers. In the latter case, we work with the mean-field disordered p spin and show that the convergence to equilibrium or to the nonergodic phase is accelerated during both the β- and α-relaxation stages. An interpretation in terms of trajectories in phase space and of an accidental fluctuation-dissipation theorem is provided.

Molecular Dynamics Simulation Studies on the Aggregation of Amyloid-β Peptides and Their Disaggregation by Ultrasonic Wave and Infrared Laser Irradiation

Alzheimer’s disease is understood to be caused by amyloid fibrils and oligomers formed by aggregated amyloid-β (Aβ) peptides. This review article presents molecular dynamics (MD) simulation studies of Aβ peptides and Aβ fragments on their aggregation, aggregation inhibition, amyloid fibril conformations in equilibrium, and disruption of the amyloid fibril by ultrasonic wave and infrared laser irradiation. In the aggregation of Aβ, a β-hairpin structure promotes the formation of intermolecular β-sheet structures. Aβ peptides tend to exist at hydrophilic/hydrophobic interfaces and form more β-hairpin structures than in bulk water. These facts are the reasons why the aggregation is accelerated at the interface. We also explain how polyphenols, which are attracting attention as aggregation inhibitors of Aβ peptides, interact with Aβ. An MD simulation study of the Aβ amyloid fibrils in equilibrium is also presented: the Aβ amyloid fibril has a different structure at one end from that at the other end. The amyloid fibrils can be destroyed by ultrasonic wave and infrared laser irradiation. The molecular mechanisms of these amyloid fibril disruptions are also explained, particularly focusing on the function of water molecules. Finally, we discuss the prospects for developing treatments for Alzheimer’s disease using MD simulations.

Replica permutation with solute tempering for molecular dynamics simulation and its application to the dimerization of amyloid-β fragments

We propose the replica permutation with solute tempering (RPST) by combining the replica-permutation method (RPM) and the replica exchange with solute tempering (REST). Temperature permutations are performed among more than two replicas in RPM, whereas temperature exchanges are performed between two replicas in the replica-exchange method (REM). The temperature transition in RPM occurs more efficiently than in REM. In REST, only the temperatures of the solute region, the solute temperatures, are exchanged to reduce the number of replicas compared to REM. Therefore, RPST is expected to be an improved method taking advantage of these methods. For comparison, we applied RPST, REST, RPM, and REM to two amyloid-β(16–22) peptides in explicit water. We calculated the transition ratio and the number of tunneling events in the temperature space and the number of dimerization events of amyloid-β(16–22) peptides. The results indicate that, in RPST, the number of replicas necessary for frequent random walks in the temperature and conformational spaces is reduced compared to the other three methods. In addition, we focused on the dimerization process of amyloid-β(16–22) peptides. The RPST simulation with a relatively small number of replicas shows that the two amyloid-β(16–22) peptides form the intermolecular antiparallel β-bridges due to the hydrophilic side-chain contact between Lys and Glu and hydrophobic side-chain contact between Leu, Val, and Phe, which stabilizes the dimer of the peptides.

All-Atom Molecular Dynamics Simulation Methods for the Aggregation of Protein and Peptides: Replica Exchange/Permutation and Nonequilibrium Simulations

Protein aggregates are associated with more than 40 serious human diseases. To understand the formation mechanism of protein aggregates at atomic level, all-atom molecular dynamics (MD) simulation is a powerful computational tool. In this chapter, we review the all-atom MD simulation methods that are useful for study on the protein aggregation. We first explain conventional MD simulation methods in physical statistical ensembles, such as the canonical and isothermal-isobaric ensembles. We then describe the generalized-ensemble algorithms such as replica-exchange and replica-permutation MD methods. These methods can overcome a difficulty, in which simulations tend to get trapped in local-minimum free-energy states. Finally we explain the nonequilibrium MD method. Some simulation results based on these methods are also presented.

Large-scale dynamics of event-chain Monte Carlo

Event-chain Monte Carlo (ECMC) accelerates the sampling of hard-sphere systems, and has been generalized to the potentials used in classical molecular simulations. Rather than imposing detailed balance on the transition probabilities, the method enforces a weaker global-balance condition in order to guarantee convergence to equilibrium. In this paper, we generalize the factor-field variant of ECMC to higher space dimensions. In the two-dimensional fluid phase, factor-field ECMC saturates the lower bound z=0 for the dynamical scaling exponent for local dynamics, whereas molecular dynamics is characterized by z=1 and local Metropolis Monte Carlo by z=2. In the presence of hexatic order, factor fields are not found to speed up the convergence. We note that generalizations of factor fields could couple to orientational order.

Diffuse optical tomography by simulated annealing via a spin Hamiltonian

Diffuse optical tomography (DOT) is an imaging modality that uses near-infrared light. Although iterative numerical schemes are commonly used for its inverse problem, correct solutions are not obtained unless good initial guesses are chosen. We propose a numerical scheme of DOT, which works even when good initial guesses of optical parameters are not available. We use simulated annealing (SA), which is a method of the Markov-chain Monte Carlo. To implement SA for DOT, a spin Hamiltonian is introduced in the cost function, and the Metropolis algorithm or single-component Metropolis-Hastings algorithm is used. By numerical experiments, it is shown that an initial random spin configuration is brought to a converged configuration by SA, and targets in the medium are reconstructed. The proposed numerical method solves the inverse problem for DOT by finding the ground state of a spin Hamiltonian with SA.

Role of current fluctuations in nonreversible samplers

It is known that the distribution of nonreversible Markov processes breaking the detailed balance condition converges faster to the stationary distribution compared to reversible processes having the same stationary distribution. This is used in practice to accelerate Markov chain Monte Carlo algorithms that sample the Gibbs distribution by adding nonreversible transitions or nongradient drift terms. The breaking of detailed balance also accelerates the convergence of empirical estimators to their ergodic expectation in the long-time limit. Here, we give a physical interpretation of this second form of acceleration in terms of currents associated with the fluctuations of empirical estimators using the level 2.5 of large deviations, which characterizes the likelihood of density and current fluctuations in Markov processes. Focusing on diffusion processes, we show that there is accelerated convergence because estimator fluctuations arise in general with current fluctuations, leading to an added large deviation cost compared to the reversible case, which shows no current. We study the current fluctuation most likely to arise in conjunction with a given estimator fluctuation and provide bounds on the acceleration, based on approximations of this current. We illustrate these results for the Ornstein-Uhlenbeck process in two dimensions and the Brownian motion on the circle.

Spatial Monte Carlo integration with annealed importance sampling

Evaluating expectations on an Ising model (or Boltzmann machine) is essential for various applications, including statistical machine learning. However, in general, the evaluation is computationally difficult because it involves intractable multiple summations or integrations; therefore, it requires approximation. Monte Carlo integration (MCI) is a well-known approximation method; a more effective MCI-like approximation method was proposed recently, called spatial Monte Carlo integration (SMCI). However, the estimations obtained using SMCI (and MCI) exhibit a low accuracy in Ising models under a low temperature owing to degradation of the sampling quality. Annealed importance sampling (AIS) is a type of importance sampling based on Markov chain Monte Carlo methods that can suppress performance degradation in low-temperature regions with the force of importance weights. In this study, a method is proposed to evaluate the expectations on Ising models combining AIS and SMCI. The proposed method performs efficiently in both high- and low-temperature regions, which is demonstrated theoretically and numerically.

Dimerization of α-Synuclein Fragments Studied by Isothermal-Isobaric Replica-Permutation Molecular Dynamics Simulation

Aggregates and fibrils of intrinsically disordered α-synuclein are associated with Parkinson's disease. Within a non-amyloid β component (NAC) spanning from the 61st to the 95th residue of α-synuclein, an 11-residue segment called NACore (₆₈GAVVTGVTAVA₇₈) is an essential region for both fibril formation and cytotoxicity. Although NACore peptides alone are known to form aggregates and amyloid fibrils, the mechanisms of aggregation and fibrillation remain unknown. This study investigated the dimerization process of NACore peptides as the initial stage of the aggregation and fibrillation processes. We performed an isothermal-isobaric replica-permutation molecular dynamics simulation, which is one of the efficient sampling methods, for the two NACore peptides in explicit water over 96 μs. The simulation succeeded in sampling a variety of dimer structures. An analysis of secondary structure revealed that most of the NACore dimers form intermolecular β-bridges. In particular, more antiparallel β-bridges were observed than parallel β-bridges. We also found that intramolecular secondary structures such as α-helix and antiparallel β-bridge are stabilized in the pre-dimer state. However, we identified that the intermolecular β-bridges tend to form directly between residues with no specific structure rather than via the intramolecular β-bridges. This is because the NACore peptides still have a low propensity to form the intramolecular secondary structures even though they are stabilized in the pre-dimer state.

Geometric allocation approach to accelerating directed worm algorithm

The worm algorithm is a versatile technique in the Markov chain Monte Carlo method for both classical and quantum systems. The algorithm substantially alleviates critical slowing down and reduces the dynamic critical exponents of various classical systems. It is crucial to improve the algorithm and push the boundary of the Monte Carlo method for physical systems. We here propose a directed worm algorithm that significantly improves computational efficiency. We use the geometric allocation approach to optimize the worm scattering process: worm backscattering is averted, and forward scattering is favored. Our approach successfully enhances the diffusivity of the worm head (kink), which is evident in the probability distribution of the relative position of the two kinks. Performance improvement is demonstrated for the Ising model at the critical temperature by measurement of exponential autocorrelation times and asymptotic variances. The present worm update is approximately 25 times as efficient as the conventional worm update for the simple cubic lattice model. Surprisingly, our algorithm is even more efficient than the Wolff cluster algorithm, which is one of the best update algorithms. We estimate the dynamic critical exponent of the simple cubic lattice Ising model to be z≈0.27 in the worm update. The worm and the Wolff algorithms produce different exponents of the integrated autocorrelation time of the magnetic susceptibility estimator but the same exponent of the asymptotic variance. We also discuss how to quantify the computational efficiency of the Markov chain Monte Carlo method. Our approach can be applied to a wide range of physical systems, such as the |ϕ|^{4} model, the Potts model, the O(n) loop model, and lattice QCD.

Promotion and Inhibition of Amyloid-β Peptide Aggregation: Molecular Dynamics Studies

Aggregates of amyloid-β (Aβ) peptides are known to be related to Alzheimer’s disease. Their aggregation is enhanced at hydrophilic–hydrophobic interfaces, such as a cell membrane surface and air-water interface, and is inhibited by polyphenols, such as myricetin and rosmarinic acid. We review molecular dynamics (MD) simulation approaches of a full-length Aβ peptide, Aβ40, and Aβ(16–22) fragments in these environments. Since these peptides have both hydrophilic and hydrophobic amino acid residues, they tend to exist at the interfaces. The high concentration of the peptides accelerates the aggregation there. In addition, Aβ40 forms a β-hairpin structure, and this structure accelerates the aggregation. We also describe the inhibition mechanism of the Aβ(16–22) aggregation by polyphenols. The aggregation of Aβ(16–22) fragments is caused mainly by the electrostatic attraction between charged amino acid residues known as Lys16 and Glu22. Since polyphenols form hydrogen bonds between their hydroxy and carboxyl groups and these charged amino acid residues, they inhibit the aggregation.

Nonequilibrium phase transition in an Ising model without detailed balance

We study a two-dimensional ferromagnetic Ising model in which spins are updated using modified versions of the Metropolis and Glauber algorithms. These update rules do not obey the detailed balance condition. The steady-state behavior of the model is studied using molecular field theory and Monte Carlo simulations. This model is found to exhibit a nonequilibrium phase transition from a "paramagnetic" state with zero magnetization to a "ferromagnetic" state with nonzero magnetization as the variable that plays the role of temperature in the spin updates is decreased. From detailed Monte Carlo simulations using the modified Metropolis algorithm, we demonstrate explicitly the nonequilibrium nature of the transition and show that it cannot be described as an equilibrium transition with an effective temperature different from the temperature used in the spin updates. The critical exponents that characterize the singular behavior near this continuous phase transition are calculated from finite size scaling of specific heat, magnetization, susceptibility, and correlation length. We find that the values of these exponents are the same (within error bars) as those of the equilibrium Ising model in two dimensions.

Simulated tempering with irreversible Gibbs sampling techniques

We present here two novel algorithms for simulated tempering simulations, which break the detailed balance condition (DBC) but satisfy the skewed detailed balance to ensure invariance of the target distribution. The irreversible methods we present here are based on Gibbs sampling and concern breaking DBC at the update scheme of the temperature swaps. We utilize three systems as a test bed for our methods: a Markov chain Monte Carlo simulation on a simple system described by a one-dimensional double well potential, the Ising model, and molecular dynamics simulations on alanine pentapeptide (ALA5). The relaxation times of inverse temperature, magnetic susceptibility, and energy density for the Ising model indicate clear gains in sampling efficiency over conventional Gibbs sampling techniques with DBC and also over the conventionally used simulated tempering with the Metropolis-Hastings (MH) scheme. Simulations on ALA5 with a large number of temperatures indicate distinct gains in mixing times for inverse temperature and consequently the energy of the system compared to conventional MH. With no additional computational overhead, our methods were found to be more efficient alternatives to the conventionally used simulated tempering methods with DBC. Our algorithms should be particularly advantageous in simulations of large systems with many temperature ladders, as our algorithms showed a more favorable constant scaling in Ising spin systems as compared with both reversible and irreversible MH algorithms. In future applications, our irreversible methods can also be easily tailored to utilize a given dynamical variable other than temperature to flatten rugged free energy landscapes.

Energetics and kinetics of substrate analog-coupled staphylococcal nuclease folding revealed by a statistical mechanical approach

Protein conformational changes associated with ligand binding, especially those involving intrinsically disordered proteins, are mediated by tightly coupled intra- and intermolecular events. Such reactions are often discussed in terms of two limiting kinetic mechanisms, conformational selection (CS), where folding precedes binding, and induced fit (IF), where binding precedes folding. It has been shown that coupled folding/binding reactions can proceed along both CS and IF pathways with the flux ratio depending on conditions such as ligand concentration. However, the structural and energetic basis of such complex reactions remains poorly understood. Therefore, we used experimental, theoretical, and computational approaches to explore structural and energetic aspects of the coupled-folding/binding reaction of staphylococcal nuclease in the presence of the substrate analog adenosine-3',5'-diphosphate. Optically monitored equilibrium and kinetic data, combined with a statistical mechanical model, gave deeper insight into the relative importance of specific and Coulombic protein-ligand interactions in governing the reaction mechanism. We also investigated structural aspects of the reaction at the residue level using NMR and all-atom replica-permutation molecular dynamics simulations. Both approaches yielded clear evidence for accumulation of a transient protein-ligand encounter complex early in the reaction under IF-dominant conditions. Quantitative analysis of the equilibrium/kinetic folding revealed that the ligand-dependent CS-to-IF shift resulted from stabilization of the compact transition state primarily by weakly ligand-dependent Coulombic interactions with smaller contributions from specific binding energies. At a more macroscopic level, the CS-to-IF shift was represented as a displacement of the reaction "route" on the free energy surface, which was consistent with a flux analysis.

Efficient Irreversible Monte Carlo Samplers

We present here two irreversible Markov chain Monte Carlo algorithms for general discrete state systems. One of the algorithms is based on the random-scan Gibbs sampler for discrete states and the other on its improved version, the Metropolized-Gibbs sampler. The algorithms we present incorporate the lifting framework with skewed detailed balance condition and construct irreversible Markov chains that satisfy the balance condition. We have applied our algorithms to 1D 4-state Potts model. The integrated autocorrelation times for magnetization and energy density indicate a reduction of the dynamical scaling exponent from z ≈ 1 to z ≈ 1/2. In addition, we have generalized an irreversible Metropolis-Hastings algorithm with skewed detailed balance, initially introduced by Turitsyn et al. [ Physica D 2011, 240, 410] for the mean field Ising model, to be now readily applicable to classical spin systems in general; application to 1D 4-state Potts model indicate a square root reduction of the mixing time at high temperatures.

Replica sub-permutation method for molecular dynamics and monte carlo simulations

We propose an improvement of the replica-exchange and replica-permutation methods, which we call the replica sub-permutation method (RSPM). Instead of considering all permutations, this method uses a new algorithm referred to as sub-permutation to perform parameter transition. The RSPM succeeds in reducing the number of combinations between replicas and parameters without the loss of sampling efficiency. For comparison, we applied the replica sub-permutation, replica-permutation, and replica-exchange methods to a β-hairpin mini protein, chignolin, in explicit water. We calculated the transition ratio and number of tunneling events in the parameter space, the number of folding-unfolding events, the autocorrelation function, and the autocorrelation time as measures of sampling efficiency. The results indicate that among the three methods, the proposed RSPM is the most efficient in both parameter and conformational spaces. © 2019 Wiley Periodicals, Inc.

Generation of ice states through deep reinforcement learning

We present a deep reinforcement learning framework where a machine agent is trained to search for a policy to generate a ground state for the square ice model by exploring the physical environment. After training, the agent is capable of proposing a sequence of local moves to achieve the goal. Analysis of the trained policy and the state value function indicates that the ice rule and loop-closing condition are learned without prior knowledge. We test the trained policy as a sampler in the Markov chain Monte Carlo and benchmark against the baseline loop algorithm. This framework can be generalized to other models with topological constraints where generation of constraint-preserving states is difficult.

Molecular simulations by generalized-ensemble algorithms in isothermal-isobaric ensemble

Generalized-ensemble algorithms are powerful techniques for investigating biomolecules such as protein, DNA, lipid membrane, and glycan. The generalized-ensemble algorithms were originally developed in the canonical ensemble. On the other hand, not only temperature but also pressure is controlled in experiments. Additionally, pressure is used as perturbation to study relationship between function and structure of biomolecules. For this reason, it is important to perform efficient conformation sampling based on the isothermal-isobaric ensemble. In this article, we review a series of the generalized-ensemble algorithms in the isothermal-isobaric ensemble: multibaric-multithermal, pressure- and temperature-simulated tempering, replica-exchange, and replica-permutation methods. These methods achieve more efficient simulation than the conventional isothermal-isobaric simulation. Furthermore, the isothermal-isobaric generalized-ensemble simulation samples conformations of biomolecules from wider range of temperature and pressure. Thus, we can estimate physical quantities more accurately at any temperature and pressure values. The applications to the biomolecular system are also presented.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Direct coupling of first-principles calculations with replica exchange Monte Carlo sampling of ion disorder in solids

We demonstrate the feasibility of performing sufficient configurational sampling of disordered oxides directly from first-principles without resorting to the use of fitted models such as cluster expansion. This is achieved by harnessing the power of modern-day cluster supercomputers using the replica exchange Monte Carlo method coupled directly with structural relaxation and energy calculation performed by density functional codes. The idea is applied successfully to the calculation of the temperature-dependence of the degree of inversion in the cation sublattice of MgAl₂O₄ spinel oxide. The possibility of bypassing fitting models will lead to investigation of disordered systems where cluster expansion is known to perform badly, for example, systems with large lattice deformation due to defects, or systems where long-range interactions dominate such as electrochemical interfaces.

Replica-Exchange Methods for Biomolecular Simulations

In this study, a replica-exchange method was developed to overcome conformational sampling difficulties in computer simulations of spin glass or other systems with rugged free-energy landscapes. This method was then applied to the protein-folding problem in combination with molecular dynamics (MD) simulation. Owing to its simplicity and sampling efficiency, the replica-exchange method has been applied to many other biological problems and has been continuously improved. The method has often been combined with other sampling techniques, such as umbrella sampling, free-energy perturbation, metadynamics, and Gaussian accelerated MD (GaMD). In this chapter, we first summarize the original replica-exchange molecular dynamics (REMD) method and discuss how new algorithms related to the original method are implemented to add new features. Heterogeneous and flexible structures of an N-glycan in a solution are simulated as an example of applications by REMD, replica exchange with solute tempering, and GaMD. The sampling efficiency of these methods on the N-glycan system and the convergence of the free-energy changes are compared. REMD simulation protocols and trajectory analysis using the GENESIS software are provided to facilitate the practical use of advanced simulation methods.

Dynamical Signature of Symmetry Fractionalization in Frustrated Magnets

The nontriviality of quantum spin liquids (QSLs) typically manifests in the nonlocal observables that signify their existence; however, this fact actually casts a shadow on detecting QSLs with experimentally accessible probes. Here, we provide a solution by unbiasedly demonstrating a dynamical signature of anyonic excitations and symmetry fractionalization in QSLs. Employing large-scale quantum Monte Carlo simulation and stochastic analytic continuation, we investigate the extended XXZ model on the kagome lattice, and find out that, across the phase transitions from Z_{2} QSLs to different symmetry breaking phases, spin spectral functions can reveal the presence and condensation of emergent anyonic spinon and vison excitations, in particular, the translational symmetry fractionalization of the latter, which can be served as the dynamical signature of the seemingly ephemeral QSLs in spectroscopic techniques such as inelastic neutron or resonance (inelastic) x-ray scatterings.

Quantum Spin Liquid with Even Ising Gauge Field Structure on Kagome Lattice

Employing large-scale quantum Monte Carlo simulations, we study the extended XXZ model on the kagome lattice. A Z_{2} quantum spin liquid phase with effective even Ising gauge field structure emerges from the delicate balance among three symmetry-breaking phases including stripe solid, staggered solid, and ferromagnet. This Z_{2} spin liquid is stabilized by an extended interaction related to the Rokhsar-Kivelson potential in the quantum dimer model limit. The phase transitions from the staggered solid to a spin liquid or ferromagnet are found to be first order and so is the transition between the stripe solid and ferromagnet. However, the transition between a spin liquid and ferromagnet is found to be continuous and belongs to the 3D XY^{*} universality class associated with the condensation of spinons. The transition between a spin liquid and stripe solid appears to be continuous and associated with the condensation of visons.

Development of isothermal-isobaric replica-permutation method for molecular dynamics and Monte Carlo simulations and its application to reveal temperature and pressure dependence of folded, misfolded, and unfolded states of chignolin

We developed a two-dimensional replica-permutation molecular dynamics method in the isothermal-isobaric ensemble. The replica-permutation method is a better alternative to the replica-exchange method. It was originally developed in the canonical ensemble. This method employs the Suwa-Todo algorithm, instead of the Metropolis algorithm, to perform permutations of temperatures and pressures among more than two replicas so that the rejection ratio can be minimized. We showed that the isothermal-isobaric replica-permutation method performs better sampling efficiency than the isothermal-isobaric replica-exchange method and infinite swapping method. We applied this method to a β-hairpin mini protein, chignolin. In this simulation, we observed not only the folded state but also the misfolded state. We calculated the temperature and pressure dependence of the fractions on the folded, misfolded, and unfolded states. Differences in partial molar enthalpy, internal energy, entropy, partial molar volume, and heat capacity were also determined and agreed well with experimental data. We observed a new phenomenon that misfolded chignolin becomes more stable under high-pressure conditions. We also revealed this mechanism of the stability as follows: TYR2 and TRP9 side chains cover the hydrogen bonds that form a β-hairpin structure. The hydrogen bonds are protected from the water molecules that approach the protein as the pressure increases.

Lifted worm algorithm for the Ising model

We design an irreversible worm algorithm for the zero-field ferromagnetic Ising model by using the lifting technique. We study the dynamic critical behavior of an energylike observable on both the complete graph and toroidal grids, and compare our findings with reversible algorithms such as the Prokof'ev-Svistunov worm algorithm. Our results show that the lifted worm algorithm improves the dynamic exponent of the energylike observable on the complete graph and leads to a significant constant improvement on toroidal grids.

Bayesian optimization for computationally extensive probability distributions

An efficient method for finding a better maximizer of computationally extensive probability distributions is proposed on the basis of a Bayesian optimization technique. A key idea of the proposed method is to use extreme values of acquisition functions by Gaussian processes for the next training phase, which should be located near a local maximum or a global maximum of the probability distribution. Our Bayesian optimization technique is applied to the posterior distribution in the effective physical model estimation, which is a computationally extensive probability distribution. Even when the number of sampling points on the posterior distributions is fixed to be small, the Bayesian optimization provides a better maximizer of the posterior distributions in comparison to those by the random search method, the steepest descent method, or the Monte Carlo method. Furthermore, the Bayesian optimization improves the results efficiently by combining the steepest descent method and thus it is a powerful tool to search for a better maximizer of computationally extensive probability distributions.

Quantum Speed Limit is Not Quantum

The quantum speed limit (QSL), or the energy-time uncertainty relation, describes the fundamental maximum rate for quantum time evolution and has been regarded as being unique in quantum mechanics. In this study, we obtain a classical speed limit corresponding to the QSL using Hilbert space for the classical Liouville equation. Thus, classical mechanics has a fundamental speed limit, and the QSL is not a purely quantum phenomenon but a universal dynamical property of Hilbert space. Furthermore, we obtain similar speed limits for the imaginary-time Schrödinger equations such as the classical master equation.

Theoretical approaches for dynamical ordering of biomolecular systems

BACKGROUND

Living systems are characterized by the dynamic assembly and disassembly of biomolecules. The dynamical ordering mechanism of these biomolecules has been investigated both experimentally and theoretically. The main theoretical approaches include quantum mechanical (QM) calculation, all-atom (AA) modeling, and coarse-grained (CG) modeling. The selected approach depends on the size of the target system (which differs among electrons, atoms, molecules, and molecular assemblies). These hierarchal approaches can be combined with molecular dynamics (MD) simulation and/or integral equation theories for liquids, which cover all size hierarchies.

SCOPE OF REVIEW

We review the framework of quantum mechanical/molecular mechanical (QM/MM) calculations, AA MD simulations, CG modeling, and integral equation theories. Applications of these methods to the dynamical ordering of biomolecular systems are also exemplified.

MAJOR CONCLUSIONS

The QM/MM calculation enables the study of chemical reactions. The AA MD simulation, which omits the QM calculation, can follow longer time-scale phenomena. By reducing the number of degrees of freedom and the computational cost, CG modeling can follow much longer time-scale phenomena than AA modeling. Integral equation theories for liquids elucidate the liquid structure, for example, whether the liquid follows a radial distribution function.

GENERAL SIGNIFICANCE

These theoretical approaches can analyze the dynamic behaviors of biomolecular systems. They also provide useful tools for exploring the dynamic ordering systems of biomolecules, such as self-assembly. This article is part of a Special Issue entitled "Biophysical Exploration of Dynamical Ordering of Biomolecular Systems" edited by Dr. Koichi Kato.

Temperature-shuffled parallel cascade selection molecular dynamics accelerates the structural transitions of proteins

Parallel cascade selection molecular dynamics (PaCS-MD) is an enhanced conformational sampling method for searching structural transition pathways from a given reactant to a product. Recently, a temperature-aided PaCS-MD (Vinod et al., Eur. Biophys. J. 2016, 45, 463) has been proposed as its extension, in which the temperatures were introduced as additional parameters in conformational resampling, whereas the temperature is fixed in the original PaCS-MD. In the present study, temperature-shuffled PaCS-MD is proposed as a further extension of temperature-aided PaCS-MD in which the temperatures are shuffled among different replicas at the beginning of each cycle of conformational resampling. To evaluate their conformational sampling efficiencies, the original, temperature-aided, and temperature-shuffled PaCS-MD were applied to a protein-folding process of Trp-cage, and their minimum computational costs to identify the native state were addressed. Through the evaluation, it was confirmed that temperature-shuffled PaCS-MD remarkably accelerated the protein-folding process of Trp-cage compared with the other methods. © 2017 Wiley Periodicals, Inc.

Enhanced conformational sampling to visualize a free-energy landscape of protein complex formation

We introduce various, recently developed, generalized ensemble methods, which are useful to sample various molecular configurations emerging in the process of protein-protein or protein-ligand binding. The methods introduced here are those that have been or will be applied to biomolecular binding, where the biomolecules are treated as flexible molecules expressed by an all-atom model in an explicit solvent. Sampling produces an ensemble of conformations (snapshots) that are thermodynamically probable at room temperature. Then, projection of those conformations to an abstract low-dimensional space generates a free-energy landscape. As an example, we show a landscape of homo-dimer formation of an endothelin-1-like molecule computed using a generalized ensemble method. The lowest free-energy cluster at room temperature coincided precisely with the experimentally determined complex structure. Two minor clusters were also found in the landscape, which were largely different from the native complex form. Although those clusters were isolated at room temperature, with rising temperature a pathway emerged linking the lowest and second-lowest free-energy clusters, and a further temperature increment connected all the clusters. This exemplifies that the generalized ensemble method is a powerful tool for computing the free-energy landscape, by which one can discuss the thermodynamic stability of clusters and the temperature dependence of the cluster networks.

Oligomer Formation of Amyloid-β(29-42) from Its Monomers Using the Hamiltonian Replica-Permutation Molecular Dynamics Simulation

Oligomers of amyloid-β peptides (Aβ) are formed during the early stage of the amyloidogenesis process and exhibit neurotoxicity. The oligomer formation process of Aβ and even that of Aβ fragments are still poorly understood, though understanding of these processes is essential for remedying Alzheimer's disease. In order to better understand the oligomerization process of the C-terminal Aβ fragment Aβ(29-42) at the atomic level, we performed the Hamiltonian replica-permutation molecular dynamics simulation with Aβ(29-42) molecules using the explicit water solvent model. We observed that oligomers increased in size through the sequential addition of monomers to the oligomer, rather than through the assembly of small oligomers. Moreover, solvent effects played an important role in this oligomerization process.

Eigenvalue analysis of an irreversible random walk with skew detailed balance conditions

An irreversible Markov-chain Monte Carlo (MCMC) algorithm with skew detailed balance conditions originally proposed by Turitsyn et al. is extended to general discrete systems on the basis of the Metropolis-Hastings scheme. To evaluate the efficiency of our proposed method, the relaxation dynamics of the slowest mode and the asymptotic variance are studied analytically in a random walk on one dimension. It is found that the performance in irreversible MCMC methods violating the detailed balance condition is improved by appropriately choosing parameters in the algorithm.

Enhanced sampling simulations to construct free-energy landscape of protein-partner substrate interaction

Molecular dynamics (MD) simulations using all-atom and explicit solvent models provide valuable information on the detailed behavior of protein-partner substrate binding at the atomic level. As the power of computational resources increase, MD simulations are being used more widely and easily. However, it is still difficult to investigate the thermodynamic properties of protein-partner substrate binding and protein folding with conventional MD simulations. Enhanced sampling methods have been developed to sample conformations that reflect equilibrium conditions in a more efficient manner than conventional MD simulations, thereby allowing the construction of accurate free-energy landscapes. In this review, we discuss these enhanced sampling methods using a series of case-by-case examples. In particular, we review enhanced sampling methods conforming to trivial trajectory parallelization, virtual-system coupled multicanonical MD, and adaptive lambda square dynamics. These methods have been recently developed based on the existing method of multicanonical MD simulation. Their applications are reviewed with an emphasis on describing their practical implementation. In our concluding remarks we explore extensions of the enhanced sampling methods that may allow for even more efficient sampling.

Conflict between fastest relaxation of a Markov process and detailed balance condition

We consider the optimization of Markovian dynamics to pursue the fastest convergence to the stationary state. The brachistochrone method is applied to the continuous-time master equation for finite-size systems. The principle of least action leads to a brachistochrone equation for the transition-rate matrix. Three-state systems are explicitly analyzed, and we find that the solution violates the detailed balance condition. The properties of the solution are studied in detail to observe the optimality of the solution. We also discuss the counterdiabatic driving for the Markovian dynamics. The transition-rate matrix is then divided into two parts, and the state is given by an eigenstate of the first part. The second part violates the detailed balance condition and plays the role of a counterdiabatic term.

Event-chain algorithm for the Heisenberg model: Evidence for z≃1 dynamic scaling

We apply the event-chain Monte Carlo algorithm to the three-dimensional ferromagnetic Heisenberg model. The algorithm is rejection-free and also realizes an irreversible Markov chain that satisfies global balance. The autocorrelation functions of the magnetic susceptibility and the energy indicate a dynamical critical exponent z≈1 at the critical temperature, while that of the magnetization does not measure the performance of the algorithm. We show that the event-chain Monte Carlo algorithm substantially reduces the dynamical critical exponent from the conventional value of z≃2.

Mass-scaling replica-exchange molecular dynamics optimizes computational resources with simpler algorithm

We develop a novel method of replica-exchange molecular dynamics (REMD) simulation, mass-scaling REMD (MSREMD) method, which improves numerical stability of simulations. In addition, the MSREMD method can also simplify a replica-exchange routine by eliminating velocity scaling. As a pilot system, a Lennard-Jones fluid is simulated with the new method. The results suggest that the MSREMD method improves the numerical stability at high temperatures compared with the conventional REMD method. For the Nosé-Hoover thermostats, we analytically demonstrate that the MSREMD simulations can reproduce completely the same trajectories of the conventional REMD ones with shorter time steps at high temperatures. Accordingly, we can easily compare the computational costs of the REMD and MSREMD simulations. We conclude that the MSREMD method decreases the instability and optimizes the computational resources with simpler algorithm.

Simulated tempering based on global balance or detailed balance conditions: Suwa-Todo, heat bath, and Metropolis algorithms

Simulated tempering (ST) is a useful method to enhance sampling of molecular simulations. When ST is used, the Metropolis algorithm, which satisfies the detailed balance condition, is usually applied to calculate the transition probability. Recently, an alternative method that satisfies the global balance condition instead of the detailed balance condition has been proposed by Suwa and Todo. In this study, ST method with the Suwa-Todo algorithm is proposed. Molecular dynamics simulations with ST are performed with three algorithms (the Metropolis, heat bath, and Suwa-Todo algorithms) to calculate the transition probability. Among the three algorithms, the Suwa-Todo algorithm yields the highest acceptance ratio and the shortest autocorrelation time. These suggest that sampling by a ST simulation with the Suwa-Todo algorithm is most efficient. In addition, because the acceptance ratio of the Suwa-Todo algorithm is higher than that of the Metropolis algorithm, the number of temperature states can be reduced by 25% for the Suwa-Todo algorithm when compared with the Metropolis algorithm.