Replica Exchange Monte Carlo Simulations Research Articles

Cluster-Spin-Glass Magnetic Behavior and Morphology in the Coordination Polymer Alloys FeyCo1-yBTT.

We report the synthesis of a series of pseudo-1D coordination polymer (CP) materials with the formula FeyCo1-yBTT (BTT = 1,3,5-benzenetrithiolate). These materials were structurally characterized by PXRD Rietveld, EXAFS, and PDF analyses, revealing that the CP superstructure enables a continuous and isomorphous alloy between the two homometallic compounds. Lower Fe loadings exhibit emergent spin glass magnetic behavior, such as memory effects and composition-dependent spin glass response time constants ranging from 6.9 × 10-9 s to 1.8 × 10-6 s. These data are consistent with the formation of spin clusters within the lattice. The magnetic behavior in these materials was modeled via replica exchange Monte Carlo simulation, which provides a good match for the experimentally measured spin glassing and magnetic phase transitions. These findings underscore how the rigid superstructure of CP and MOF scaffolds can enable the systematic tuning of physical properties, such as the spin glass behavior described here.

Interaction models and configurational entropies of binary MoTa and the MoNbTaW high entropy alloy

We introduce a simplified method to model the interatomic interactions of high entropy alloys based on a lookup table of cluster energies. These interactions are employed in replica exchange Monte Carlo simulations with histogram analysis to obtain thermodynamic properties across a broad temperature range. Kikuchi's cluster variation method entropy formalism is applied to directly calculate entropy from statistics on short- and long-range chemical order, and we discuss the convergence of the entropy as clusters of differing size are included. A high temperature series expansion aids in our understanding of the convergence. Computer codes implementing these methods, and supporting data, are freely available on the internet.

Probing Local Environments of Oxygen Vacancies Responsible for Hydration in Sc-Doped Barium Zirconates at Elevated Temperatures: In Situ X-ray Absorption Spectroscopy, Thermogravimetry, and Active Learning Ab Initio Replica Exchange Monte Carlo Simulations

Proton-conducting oxides, specifically heavily Sc-doped barium zirconate perovskite, have attracted attention as electrolytes for intermediate-temperature protonic ceramic fuel cells because of their high proton conductivity and high chemical stability against carbon dioxide in that temperature regime. Hydration is a key reaction for incorporating protons by filling oxygen vacancies, VO, with hydroxyl groups and activating proton conduction in the perovskite. However, probing the local environment of oxygen vacancies responsible for hydration is challenging because the behavior depends on the temperature and water partial pressure, which necessitates in situ observations and calculations of the local environments at elevated temperatures. To obtain such information, we combined in situ X-ray absorption spectroscopy (XAS) for both the Sc and Zr K-edges, thermogravimetry, X-ray diffractometry, and active learning ab initio replica exchange Monte Carlo (RXMC) simulations in undoped and 20–40 at % Sc-doped barium zirconates at and below 800 °C. The presence of oxygen vacancies adjacent to Sc and Zr in the dehydrated samples and the hydration of these oxygen vacancies under a wet atmosphere were probed by in situ XAS for Sc and Zr pre-edges at elevated temperatures. Here, the microscopic hydration linearly responds to the macroscopic degree of hydration. RXMC sampling further supports the presence of Sc-VO-Zr and Sc-VO-Sc environments. An initial hydration occurs in the Sc-VO-Zr environment at and above 600 °C, but the Sc-VO-Sc environment contribution is greater at higher degrees of hydration. The Zr-Vo-Zr environment is the least abundant among them for the whole temperature range examined and thus has a negligible impact.

De novo protein fold design through sequence-independent fragment assembly simulations

De novo protein design generally consists of two steps, including structure and sequence design. Many protein design studies have focused on sequence design with scaffolds adapted from native structures in the PDB, which renders novel areas of protein structure and function space unexplored. We developed FoldDesign to create novel protein folds from specific secondary structure (SS) assignments through sequence-independent replica-exchange Monte Carlo (REMC) simulations. The method was tested on 354 non-redundant topologies, where FoldDesign consistently created stable structural folds, while recapitulating on average 87.7% of the SS elements. Meanwhile, the FoldDesign scaffolds had well-formed structures with buried residues and solvent-exposed areas closely matching their native counterparts. Despite the high fidelity to the input SS restraints and local structural characteristics of native proteins, a large portion of the designed scaffolds possessed global folds completely different from natural proteins in the PDB, highlighting the ability of FoldDesign to explore novel areas of protein fold space. Detailed data analyses revealed that the major contributions to the successful structure design lay in the optimal energy force field, which contains a balanced set of SS packing terms, and REMC simulations, which were coupled with multiple auxiliary movements to efficiently search the conformational space. Additionally, the ability to recognize and assemble uncommon super-SS geometries, rather than the unique arrangement of common SS motifs, was the key to generating novel folds. These results demonstrate a strong potential to explore both structural and functional spaces through computational design simulations that natural proteins have not reached through evolution.

Cohesin and CTCF complexes mediate contacts in chromatin loops depending on nucleosome positions

The spatial organization of the eukaryotic genome plays an important role in regulating transcriptional activity. In the nucleus, chromatin forms loops that assemble into fundamental units called topologically associating domains that facilitate or inhibit long-range contacts. These loops are formed and held together by the ring-shaped cohesin protein complex, and this can involve binding of CCCTC-binding factor (CTCF). High-resolution conformation capture experiments provide the frequency at which two DNA fragments physically associate in three-dimensional space. However, technical limitations of this approach, such as low throughput, low resolution, or noise in contact maps, make data interpretation and identification of chromatin intraloop contacts, e.g., between distal regulatory elements and their target genes, challenging. Herein, an existing coarse-grained model of chromatin at single-nucleosome resolution was extended by integrating potentials describing CTCF and cohesin. We performed replica-exchange Monte Carlo simulations with regularly spaced nucleosomes and experimentally determined nucleosome positions in the presence of cohesin-CTCF, as well as depleted systems as controls. In fully extruded loops caused by the presence of cohesin and CTCF, the number of contacts within the formed loops was increased. The number and types of these contacts were impacted by the nucleosome distribution and loop size. Microloops were observed within cohesin-mediated loops due to thermal fluctuations without additional influence of other factors, and the number, size, and shape of microloops were determined by nucleosome distribution and loop size. Nucleosome positions directly affect the spatial structure and contact probability within a loop, with presumed consequences for transcriptional activity.

The role of the second virial coefficient in the vapor-liquid phase coexistence of anisotropic square-well particles

We examine the role of the second virial coefficient in the vapor-liquid (VL) phase coexistence of anisotropic hard bodies with square-well attractions of variable range. According to the extended law of corresponding states, the parameters of hard body interactions and attractions can be built into the reduced density and the reduced second virial coefficients, respectively, which gives rise to the collapse of all VL binodals in the reduced second virial coefficient vs. reduced density plane. The second virial perturbation theory shows that the shape dependence appears as an extra parameter in the phase behavior of anisotropic particles, which does not make possible the perfect collapse of the VL binodals for varying shapes. Interestingly, the binodal curves go closely together and even cross each other in the liquid side allowing to define a quasi-master curve. The existence of an almost perfect master curve is confirmed by replica-exchange Monte Carlo simulations for oblate square-well ellipsoids with several shape anisotropies.

Optimal Replica-Exchange Molecular Simulations in Combination with Evolution Strategies.

We have incorporated Evolution Strategies into the Replica-Exchange Monte Carlo simulation method to predict the phase behavior of several example fluids. The replica-exchange method allows one system to exchange temperatures with its neighbors to search for the most stable structure relatively efficiently in a single simulation. However, if the temperature intervals of the replicas are not positioned carefully, there is an issue that local exchange does not occur. Our results for a simple Lennard-Jones fluid and the liquid-crystal Yukawa model demonstrate the utility of the approach when compared to conventional methods. When Evolution Strategies were applied to the Replica-Exchange Monte Carlo simulation, the problem of a significant localized decrease in exchange probability near the phase transition was avoided. By obtaining the optimal temperature intervals, the system efficiently traverses a broader parameter space with a small number of replicas. This is equivalent to accelerating molecular simulations with limited computational resources and can be useful when attempting to predict the phase behavior of complex systems.

GPU-I-TASSER: a GPU accelerated I-TASSER protein structure prediction tool.

Accurate and efficient predictions of protein structures play an important role in understanding their functions. Iterative Threading Assembly Refinement (I-TASSER) is one of the most successful and widely used protein structure prediction methods in the recent community-wide CASP experiments. Yet, the computational efficiency of I-TASSER is one of the limiting factors that prevent its application for large-scale structure modeling. We present I-TASSER for Graphics Processing Units (GPU-I-TASSER), a GPU accelerated I-TASSER protein structure prediction tool for fast and accurate protein structure prediction. Our implementation is based on OpenACC parallelization of the replica-exchange Monte Carlo simulations to enhance the speed of I-TASSER by extending its capabilities to the GPU architecture. On a benchmark dataset of 71 protein structures, GPU-I-TASSER achieves on average a 10× speedup with comparable structure prediction accuracy compared to the CPU version of the I-TASSER. The complete source code for GPU-I-TASSER can be downloaded and used without restriction from https://zhanggroup.org/GPU-I-TASSER/. Supplementary data are available at Bioinformatics online.

Knots are Generic Stable Phases in Semiflexible Polymers

Semiflexible polymer models are widely used as a paradigm to understand structural phases in biomolecules including folding of proteins. Since stable knots are not so common in real proteins, the existence of stable knots in semiflexible polymers has not been explored much. Here, via extensive replica exchange Monte Carlo simulation we investigate the same for a bead-stick and a bead-spring homopolymer model that covers the whole range from flexible to stiff. We establish the fact that the presence of stable knotted phases in the phase diagram is dependent on the ratio $r_b/r_{\rm{min}}$ where $r_b$ is the equilibrium bond length and $r_{\rm{min}}$ is the distance for the strongest nonbonded contacts. Our results provide evidence for both models that if the ratio $r_b/r_{\rm{min}}$ is outside a small window around unity then depending on the bending stiffness one always encounters stable knotted phases along with the usual frozen and bent-like structures at low temperatures. These findings prompt us to conclude that knots are generic stable phases in semiflexible polymers.

Anisotropy-independent packing of confined hard ellipses

We have examined the close-packed structure and ordering properties of hard ellipses between two parallel hard walls by using geometrical conjectures and Monte Carlo simulations. Starting from closely packed hard disk structures, we discover two competing densest packings of ellipses with the stretching of hard disks into different directions. We find that parallel and tilted close-packed structures alternate with widening pore width up to the two-dimensional bulk limit, where both structures produce the same close packing value (η=π/12). Our results highlight that the densest packing does not depend on the aspect ratio of the ellipse. This universal behavior is confirmed for three different aspect ratios (κ=2, 3, and 4) using replica exchange Monte Carlo simulation. Our simulation results show that even the global phase diagram is universal, where layering transitions between parallel structures and tilted-parallel structural changes are present.

Phase Transitions and Hysteresis for a Simple Model Liquid Crystal by Replica-Exchange Monte Carlo Simulations.

In this work, the advantages of applying the temperature and pressure replica-exchange method to investigate the phase transitions and the hysteresis for liquid-crystal fluids were demonstrated. In applying this method to the commonly used Hess–Su liquid-crystal model, heat capacity peaks and points of phase co-existence were observed. The absence of a smectic phase at higher densities and a narrow range of the nematic phase were reported. The identity of the crystalline phase of this system was found to a hexagonal close-packed solid. Since the nematic-solid phase transition is strongly first order, care must be taken when using this model not to inadvertently simulate meta-stable nematic states at higher densities. In further analysis, the Weighted Histogram Analysis Method was applied to verify the precise locations of the phase transition points.

Comparison of Conformational Phase Behavior for Flexible and Semiflexible Polymers.

We employ the recently introduced generalized microcanonical inflection point method for the statistical analysis of phase transitions in flexible and semiflexible polymers and study the impact of the bending stiffness upon the character and order of transitions between random-coil, globules, and pseudocrystalline conformations. The high-accuracy estimates of the microcanonical entropy and its derivatives required for this study were obtained by extensive replica-exchange Monte Carlo simulations. We observe that the transition behavior into the compact phases changes qualitatively with increasing bending stiffness. Whereas the collapse transition is less affected, the first-order liquid-solid transition characteristic for flexible polymers ceases to exist once bending effects dominate over attractive monomer-monomer interactions.

Three-step melting of hard superdisks in two dimensions.

We explore the link between the melting scenarios of two-dimensional systems of hard disks and squares through replica-exchange Monte Carlo simulations of hard superdisks. The well-known melting scenarios are observed in the disk and square limits, while we observe an unusual three-step scenario for dual shapes. We find that two mesophases mediate the melting: a hexatic phase and another fluid phase with a D_{2} local symmetry, we call it rhombatic, where both bond and particle orientational orders are quasi-long-range. Our results show that not only can the melting process of liquid-crystal forming molecules be complicated, where elongated shapes stabilize several mesophases, but also that of anisotropic quasispherical molecules.

EDock: blind protein\u2013ligand docking by replica-exchange monte carlo simulation

Protein–ligand docking is an important approach for virtual screening and protein function annotation. Although many docking methods have been developed, most require a high-resolution crystal structure of the receptor and a user-specified binding site to start. This information is, however, not available for the majority of unknown proteins, including many pharmaceutically important targets. Developing blind docking methods without predefined binding sites and working with low-resolution receptor models from protein structure prediction is thus essential. In this manuscript, we propose a novel Monte Carlo based method, EDock, for blind protein–ligand docking. For a given protein, binding sites are first predicted by sequence-profile and substructure-based comparison searches with initial ligand poses generated by graph matching. Next, replica-exchange Monte Carlo (REMC) simulations are performed for ligand conformation refinement under the guidance of a physical force field coupled with binding-site distance constraints. The method was tested on two large-scale datasets containing 535 protein–ligand pairs. Without specifying binding pockets on the experimental receptor structures, EDock achieves on average a ligand RMSD of 2.03 Å, which compares favorably with state-of-the-art docking methods including DOCK6 (2.68 Å) and AutoDock Vina (3.92 Å). When starting with predicted models from I-TASSER, EDock still generates reasonable docking models, with a success rate 159% and 67% higher than DOCK6 and AutoDock Vina, respectively. Detailed data analyses show that the major advantage of EDock lies in reliable ligand binding site predictions and extensive REMC sampling, which allows for the implementation of multiple van der Waals weightings to accommodate different levels of steric clashes and cavity distortions and therefore enhances the robustness of low-resolution docking with predicted protein structures.

Increasing the degree of parallelism using speculative execution in task-based runtime systems.

Task-based programming models have demonstrated their efficiency in the development of scientific applications on modern high-performance platforms. They allow delegation of the management of parallelization to the runtime system (RS), which is in charge of the data coherency, the scheduling, and the assignment of the work to the computational units. However, some applications have a limited degree of parallelism such that no matter how efficient the RS implementation, they may not scale on modern multicore CPUs. In this paper, we propose using speculation to unleash the parallelism when it is uncertain if some tasks will modify data, and we formalize a new methodology to enable speculative execution in a graph of tasks. This description is partially implemented in our new C++ RS called SPETABARU, which is capable of executing tasks in advance if some others are not certain to modify the data. We study the behavior of our approach to compute Monte Carlo and replica exchange Monte Carlo simulations.

EvoDesign: Designing Protein–Protein Binding Interactions Using Evolutionary Interface Profiles in Conjunction with an Optimized Physical Energy Function

EvoDesign (https://zhanglab.ccmb.med.umich.edu/EvoDesign) is an online server system for protein design. The method uses evolutionary profiles to guide the sequence search simulation and demonstrated significant advantages over physics-based approaches in terms of more accurately designing proteins that adopt desired target folds. Despite the success, the previous EvoDesign program focused only on monomer protein design, which limited its ability and usefulness in terms of designing functional proteins. In this work, we propose a new EvoDesign server, which extends the principles of evolution-based design to design protein–protein interactions. Starting from a two-chain complex structure, structurally similar interfaces are identified from known protein–protein interaction databases. An interface evolutionary profile is then constructed from a multiple sequence alignment of the interface analogies, which is combined with a newly developed, atomic-level physical energy function to guide the replica-exchange Monte Carlo simulation search. The purpose of the server is to redesign the specified complex chain to increase its stability and binding affinity for the other chain in the complex. With the improved scope and accuracy of the methodology, the new EvoDesign pipeline should become a useful online tool for functional protein design and drug discovery studies.

Influence of bonded interactions on structural phases of flexible polymers.

We introduce a novel coarse-grained bead-spring model for flexible polymers to systematically examine the effects of an adjusted bonded potential on the formation and stability of structural macrostates in a thermal environment. The density of states obtained in advanced replica-exchange Monte Carlo simulations is analyzed by employing the recently developed generalized microcanonical inflection-point analysis method, which enables the identification of diverse structural phases and the construction of a suitably parameterized hyperphase diagram. It reveals that icosahedral phases dominate for polymers with asymmetric and narrow bond potentials, whereas polymers with symmetric and more elastic bonds tend to form amorphous structures with non-icosahedral cores. We also observe a hierarchy in the freezing transition behavior associated with the formation of the surface layer after nucleation.

Modeling the Role of Excluded Volume in Zeolite Structure Direction.

We investigate the formation of zeolite structures in replica-exchange Monte Carlo simulations of a reactive model of silica polymerization. The simulations incorporate hard spheres to model the effect of excluded volume caused by structure-directing agents (SDAs). We focus on modeling the formation of cage-type zeolite frameworks SOD and LTA. Our model predicts that a relatively wide range of SDA sizes could be used to construct SOD, whereas a narrower range will work for constructing LTA. We also predict that there is potential benefit of including multiple SDAs in each zeolite unit cell, and in the case of LTA with both small and large cavities, there is a strong potential benefit using both small and large SDAs that match the cavities' sizes. We hypothesize that the volume exclusion reduces the configuration space available to the assembling silica units, making it easier for the system to find ordered structures with quasi-spherical cavities.

MR-REX: molecular replacement by cooperative conformational search and occupancy optimization on low-accuracy protein models.

Molecular replacement (MR) has commonly been employed to derive the phase information in protein crystal X-ray diffraction, but its success rate decreases rapidly when the search model is dissimilar to the target. MR-REX has been developed to perform an MR search by replica-exchange Monte Carlo simulations, which enables cooperative rotation and translation searches and simultaneous clash and occupancy optimization. MR-REX was tested on a set of 1303 protein structures of different accuracies and successfully placed 699 structures at positions that have an r.m.s.d. of below 2 Å to the target position, which is 10% higher than was obtained by Phaser. However, cases studies show that many of the models for which Phaser failed and MR-REX succeeded can be solved by Phaser by pruning them and using nondefault parameters. The factors effecting success and the parts of the methodology which lead to success are studied. The results demonstrate a new avenue for molecular replacement which outperforms (and has results that are complementary to) the state-of-the-art MR methods, in particular for distantly homologous proteins.

A Probabilistic Framework for Constructing Temporal Relations in Replica Exchange Molecular Trajectories.

Knowledge of the structure and dynamics of biomolecules is essential for elucidating the underlying mechanisms of biological processes. Given the stochastic nature of many biological processes, like protein unfolding, it is almost impossible that two independent simulations will generate the exact same sequence of events, which makes direct analysis of simulations difficult. Statistical models like Markov chains, transition networks, etc. help in shedding some light on the mechanistic nature of such processes by predicting long-time dynamics of these systems from short simulations. However, such methods fall short in analyzing trajectories with partial or no temporal information, for example, replica exchange molecular dynamics or Monte Carlo simulations. In this work, we propose a probabilistic algorithm, borrowing concepts from graph theory and machine learning, to extract reactive pathways from molecular trajectories in the absence of temporal data. A suitable vector representation was chosen to represent each frame in the macromolecular trajectory (as a series of interaction and conformational energies), and dimensionality reduction was performed using principal component analysis (PCA). The trajectory was then clustered using a density-based clustering algorithm, where each cluster represents a metastable state on the potential energy surface (PES) of the biomolecule under study. A graph was created with these clusters as nodes with the edges learned using an iterative expectation maximization algorithm. The most reactive path is conceived as the widest path along this graph. We have tested our method on RNA hairpin unfolding trajectory in aqueous urea solution. Our method makes the understanding of the mechanism of unfolding in the RNA hairpin molecule more tractable. As this method does not rely on temporal data, it can be used to analyze trajectories from Monte Carlo sampling techniques and replica exchange molecular dynamics (REMD).