Ergodic convergence properties of supercooled liquids and glasses

Quantifying the energy landscape in weakly and strongly disordered frictional media

We investigate the "roughness" of the energy landscape of a system that diffuses in a heterogeneous medium with a random position-dependent friction coefficient α(x). This random friction acting on the system stems from spatial inhomogeneity in the surrounding medium and is modeled using the generalized Caldira-Leggett model. For a weakly disordered medium exhibiting a Gaussian random diffusivity D(x) = kBT/α(x) characterized by its average value ⟨D(x)⟩ and a pair-correlation function ⟨D(x1)D(x2)⟩, we find that the renormalized intrinsic diffusion coefficient is lower than the average one due to the fluctuations in diffusivity. The induced weak internal friction leads to increased roughness in the energy landscape. When applying this idea to diffusive motion in liquid water, the dissociation energy for a hydrogen bond gradually approaches experimental findings as fluctuation parameters increase. Conversely, for a strongly disordered medium (i.e., ultrafast-folding proteins), the energy landscape ranges from a few to a few kcal/mol, depending on the strength of the disorder. By fitting protein folding dynamics to the escape process from a metastable potential, the decreased escape rate conceptualizes the role of strong internal friction. Studying the energy landscape in complex systems is helpful because it has implications for the dynamics of biological, soft, and active matter systems.

Ergodic Measure and Potential Control of Anomalous Diffusion

In statistical mechanics, the ergodic hypothesis (i.e., the long-time average is the same as the ensemble average) accompanying anomalous diffusion has become a continuous topic of research, being closely related to irreversibility and increasing entropy. While measurement time is finite for a given process, the time average of an observable quantity might be a random variable, whose distribution width narrows with time, and one wonders how long it takes for the convergence rate to become a constant. This is also the premise of ergodic establishment, because the ensemble average is always equal to the constant. We focus on the time-dependent fluctuation width for the time average of both the velocity and kinetic energy of a force-free particle described by the generalized Langevin equation, where the stationary velocity autocorrelation function is considered. Subsequently, the shortest time scale can be estimated for a system transferring from a stationary state to an effective ergodic state. Moreover, a logarithmic spatial potential is used to modulate the processes associated with free ballistic diffusion and the control of diffusion, as well as the minimal realization of the whole power-law regime. The results presented suggest that non-ergodicity mimics the sparseness of the medium and reveals the unique role of logarithmic potential in modulating diffusion behavior.

Finite-Size Effects in Simulations of Peptide/Lipid Assembly

Molecular dynamics simulations are an attractive tool for understanding lipid/peptide self-assembly but can be plagued by inaccuracies when the system sizes are too small. The general guidance from self-assembly simulations of homogeneous micelles is that the total number of surfactants should be three to five times greater than the equilibrium aggregate number of surfactants per micelle. Herein, the heuristic is tested on the more complicated self-assembly of lipids and amphipathic peptides using the Cooke and Martini 3 coarse-grained models. Cooke model simulations with 50 to 1000 lipids and no peptide are dominated by finite-size effects, with usually one aggregate (micelle or nanodisc) containing most of the lipids forming at each system size. Approximately 200 systems of different peptide/lipid (P/L) ratios and sizes of up to 1000 lipids yield a "finite-size phase diagram" for peptide driven self-assembly, including a coexistence region of micelles and discs. Insights from the Cooke model are applied to the assembly of dimyristoylphosphatidylcholine and the ELK-neutral peptide using the Martini 3 model. Systems of 150, 450, and 900 lipids with P/L = 1/6.25 form mixtures of lipid-rich discs that agree in size with experiment and peptide-rich micelles. Only the 150-lipid system shows finite-size effects, which arise from the long-tailed distribution of aggregate sizes. The general rule of three to five times the equilibrium aggregate size remains a practical heuristic for the Cooke and Martini 3 systems investigated here.

Inertia triggers nonergodicity of fractional Brownian motion

How related are the ergodic properties of the over- and underdamped Langevin equations driven by fractional Gaussian noise? We here find that for massive particles performing fractional Brownian motion (FBM) inertial effects not only destroy the stylized fact of the equivalence of the ensemble-averaged mean-squared displacement (MSD) to the time-averaged MSD (TAMSD) of overdamped or massless FBM, but also dramatically alter the values of the ergodicity-breaking parameter (EB). Our theoretical results for the behavior of EB for underdamped or massive FBM for varying particle mass m, Hurst exponent H, and trace length T are in excellent agreement with the findings of stochastic computer simulations. The current results can be of interest for the experimental community employing various single-particle-tracking techniques and aiming at assessing the degree of nonergodicity for the recorded time series (studying, e.g., the behavior of EB versus lag time). To infer FBM as a realizable model of anomalous diffusion for a set single-particle-tracking data when massive particles are being tracked, the EBs from the data should be compared to EBs of massive (rather than massless) FBM.

Simulation of a generalized asset exchange model with economic growth and wealth distribution

The agent-based yard-sale model of wealth inequality is generalized to incorporate exponential economic growth and its distribution. The distribution of economic growth is nonuniform and is determined by the wealth of each agent and a parameter λ. Our numerical results indicate that the model has a critical point at λ=1 between a phase for λ<1 with economic mobility and exponentially growing wealth of all agents and a nonstationary phase for λ≥1 with wealth condensation and no mobility. We define the energy of the system and show that the system can be considered to be in thermodynamic equilibrium for λ<1. Our estimates of various critical exponents are consistent with a mean-field theory [see W. Klein et al., following paper, Phys. Rev. E 104, 014151 (2021)10.1103/PhysRevE.104.014151]. The exponents do not obey the usual scaling laws unless a combination of parameters that we refer to as the Ginzburg parameter is held fixed as the phase transition is approached. The model illustrates that both poorer and richer agents benefit from economic growth if its distribution does not favor the richer agents too strongly. This work and the following theoretical paper contribute to our understanding of whether the methods of equilibrium statistical mechanics can be applied to economic systems.

Finite-Size Effects and Optimal System Sizes in Simulations of Surfactant Micelle Self-Assembly

The spontaneous formation of micelles in aqueous solutions is governed by the amphipathic nature of surfactants and is practically interesting due to the regular use of micelles as membrane mimics, for the characterization of protein structure, and for drug design and delivery. We performed a systematic characterization of the finite-size effect observed in single-component dodecylphosphocholine (DPC) micelles with the coarse-grained MARTINI model. Of multiple coarse-grained solvent models investigated using large system sizes, the nonpolarizable solvent model was found to most accurately reproduce SANS spectra of 100 mM DPC in aqueous solution. We systematically investigated the finite-size effect at constant 100 mM concentration in 23 systems of sizes 40-150 DPC, confirming the finite-size effect to manifest as an oscillation in the mean micelle aggregation number about the thermodynamic aggregation number as the system size increases. The oscillations in aggregation number mostly diminish once the system supports the formation of three micelles. Similar oscillations were observed in the estimated critical micelle concentration with a mean value of 1.10 mM, which is in agreement with experiment to 0.1 mM. The accuracy of using a multiscale simulation approach to avoid finite-size effects in the micelle size distribution and SANS spectra using MARTINI and CHARMM36 was explored using multiple long time scale 500 DPC coarse-grained simulations, which were back-mapped to CHARMM36 all-atom systems. It was found that the MARTINI model generally occupies more volume than the all-atom model, leading to the formation of micelles that are of a reasonable radius of gyration but are smaller in aggregation number. The systematic characterization of the finite-size effect and exploration of multiscale modeling presented in this work provide guidance for the accurate modeling of micelles in simulations.

Effective Diffusivity for Transport with Fluctuating Drift Velocity

Consider a particle whose drift velocity fluctuates due to transitions among discrete states or due to diffusion in a confined moving fluid. At long times, the dynamics of the particle in the direction of transport can be described in terms of the average drift velocity and an effective diffusivity. For both types of fluctuations, we show that the effective diffusivity is the sum of the average intrinsic diffusivity and the time integral of the velocity correlation function of the deviation of the fluctuating velocity from its mean value. For nearest-neighbor interstate transitions and for one-dimensional diffusion in a perpendicular direction, the time integral can be found in closed form. Our analytical expressions for the effective diffusivity recover the classic results for Taylor dispersion in the laminar flow of viscous fluid and for cargo transport along microtubules by molecular motors.

Parallel Prefetching for Canonical Ensemble Monte Carlo Simulations

In order to enable large-scale molecular simulations, algorithms must efficiently utilize multicore processors that continue to increase in total core count over time with relatively stagnant clock speeds. Although parallelized molecular dynamics (MD) software has taken advantage of this trend in computer hardware, single-particle perturbations with Monte Carlo (MC) are more difficult to parallelize than system-wide updates in MD using domain decomposition. Instead, prefetching reconstructs the serial Markov chain after computing multiple MC trials in parallel. Canonical ensemble MC simulations of a Lennard-Jones fluid with prefetching resulted in up to a factor of 1.7 speedup using 2 threads, and a factor of 3 speedup using 4 threads. Strategies for maximizing efficiency of prefetching simulations are discussed, including the potentially counterintuitive benefit of reduced acceptance probabilities. Determination of the optimal acceptance probability for a parallel simulation is simplified by theoretical prediction from serial simulation data. Finally, complete open-source code for parallel prefetch simulations was made available in the Free Energy and Advance Sampling Simulation Toolkit (FEASST).

Effective ergodicity breaking phase transition in a driven-dissipative system

We show that the Olami-Feder-Christensen model exhibits an effective ergodicity breaking transition as the noise is varied. Above the critical noise, the system is effectively ergodic because the time-averaged stress on each site converges to the global spatial average. In contrast, below the critical noise, the stress on individual sites becomes trapped in different limit cycles, and the system is not ergodic. To characterize this transition, we use ideas from the study of dynamical systems and compute recurrence plots and the recurrence rate. The order parameter is identified as the recurrence rate averaged over all sites and exhibits a jump at the critical noise. We also use ideas from percolation theory and analyze the clusters of failed sites to find numerical evidence that the transition, when approached from above, can be characterized by exponents that are consistent with hyperscaling.

Diffusion Monte Carlo simulations of gas phase and adsorbed D2-(H2)n clusters

We have computed ground state energies and analyzed radial distributions for several gas phase and adsorbed D₂(H₂)_n and HD(H₂)_n clusters. An external model potential designed to mimic ionic adsorption sites inside porous materials is used [M. Mella and E. Curotto, J. Phys. Chem. A 121, 5005 (2017)]. The isotopic substitution lowers the ground state energies by the expected amount based on the mass differences when these are compared with the energies of the pure clusters in the gas phase. A similar impact is found for adsorbed aggregates. The dissociation energy of D₂ from the adsorbed clusters is always much higher than that of H₂ from both pure and doped aggregates. Radial distributions of D₂ and H₂ are compared for both the gas phase and adsorbed species. For the gas phase clusters, two types of hydrogen-hydrogen interactions are considered: one based on the assumption that rotations and translations are adiabatically decoupled and the other based on nonisotropic four-dimensional potential. In the gas phase clusters of sufficiently large size, we find the heavier isotopomer more likely to be near the center of mass. However, there is a considerable overlap among the radial distributions of the two species. For the adsorbed clusters, we invariably find the heavy isotope located closer to the attractive interaction source than H₂, and at the periphery of the aggregate, H₂ molecules being substantially excluded from the interaction with the source. This finding rationalizes the dissociation energy results. For D₂-(H₂)_n clusters with n≥12, such preference leads to the desorption of D₂ from the aggregate, a phenomenon driven by the minimization of the total energy that can be obtained by reducing the confinement of (H₂)₁₂. The same happens for (H₂)₁₃, indicating that such an effect may be quite general and impact on the absorption of quantum species inside porous materials.

Structural dynamics and inhibitor searching for Wnt-4 protein using comparative computational studies

Wnt-4 (wingless mouse mammary tumor virus integration site-4) protein is involved in many crucial embryonic pathways regulating essential processes. Aberrant Wnt-4 activity causes various anomalies leading to gastric, colon, or breast cancer. Wnt-4 is a conserved protein in structure and sequence. All Wnt proteins contain an unusual fold comprising of a thumb (or N-terminal domain) and index finger (or C-terminal domain) bifurcated by a palm domain. The aim of this study was to identify the best inhibitors of Wnt-4 that not only interact with Wnt-4 protein but also with the covalently bound acyl group to inhibit aberrant Wnt-4 activity. A systematic computational approach was used to analyze inhibition of Wnt-4. Palmitoleic acid was docked into Wnt-4 protein, followed by ligand-based virtual screening of nearly 209,847 compounds; conformer generation of 271 compounds resulted from extensive virtual screening and comparative docking of 10,531 conformers of 271 unique compounds through GOLD (Genetic Optimization for Ligand Docking), AutoDock-Vina, and FRED (Fast Rigid Exhaustive Docking) was subsequently performed. Linux scripts was used to handle the libraries of compounds. The best compounds were selected on the basis of having maximum interactions to protein with bound palmitoleic acid. These represented lead inhibitors in further experiments. Palmitoleic acid is important for efficient Wnt activity, but aberrant Wnt-4 expression can be inhibited by designing inhibitors interacting with both protein and palmitoleic acid.

Hidden complexity in the isomerization dynamics of Holliday junctions

A plausible consequence of the rugged folding energy landscapes inherent to biomolecules is that there may be more than one functionally competent folded state. Indeed, molecule-to-molecule variations in the folding dynamics of enzymes and ribozymes have recently been identified in single-molecule experiments, but without systematic quantification or an understanding of their structural origin. Here, using concepts from glass physics and complementary clustering analysis, we provide a quantitative method to analyse single-molecule fluorescence resonance energy transfer (smFRET) data, thereby probing the isomerization dynamics of Holliday junctions, which display such heterogeneous dynamics over a long observation time (T(obs) ≈ 40 s). We show that the ergodicity of Holliday junction dynamics is effectively broken and that their conformational space is partitioned into a folding network of kinetically disconnected clusters. Theory suggests that the persistent heterogeneity of Holliday junction dynamics is a consequence of internal multiloops with varying sizes and flexibilities frozen by Mg(2+) ions. An annealing experiment using Mg(2+) pulses lends support to this idea by explicitly showing that interconversions between trajectories with different patterns can be induced.

Balancing simulation accuracy and efficiency with the Amber united atom force field

We have analyzed the quality of a recently proposed Amber united-atom model and its overall efficiency in ab initio folding and thermodynamic sampling of two stable beta-hairpins. It is found that the mean backbone structures are quite consistent between the simulations in the united-atom and its corresponding all-atom models in Amber. More importantly, the simulated beta turns are also consistent between the two models. Finally, the chemical shifts on H alpha are highly consistent between simulations in the two models, although the simulated chemical shifts are lower than experiment, indicating less structured peptides, probably due to the omission of the hydrophobic term in the simulations. More interestingly, the stabilities of both beta-hairpins at room temperature are similar to those derived from the NMR measurement, whether the united-atom or the all-atom model is used. Detailed analysis shows high percentages of backbone torsion angles within the beta region and high percentages of native contacts. Given the reasonable quality of the united-atom model with respect to experimental data, we have further studied the simulation efficiency of the united-atom model over the all-atom model. Our data shows that the united-atom model is a factor of 6-8 faster than the all-atom model as measured with the ab initio first pass folding time for the two tested beta-hairpins. Detailed structural analysis shows that all ab initio folded trajectories enter the native basin, whether the united-atom model or the all-atom model is used. Finally, we have also studied the simulation efficiency of the united-atom model as measured in terms of how fast thermodynamic convergence can be achieved. It is apparent that the united-atom simulations reach convergence faster than the all-atom simulations with respect to both mean potential energies and mean native contacts. These findings show that the efficiency of the united-atom model is clearly beyond the per-step dynamics simulation of about 2 over the all-atom model. Thus, reasonable reduction of a protein model can be achieved with improved sampling efficiency while still preserving a high level of accuracy for applications in both ab initio folding and thermodynamic sampling. This study motivates us to develop more simplified protein models with sufficient consistency with the all-atom models for enhanced conformational sampling.

Cellular automaton model of damage

We investigate the role of equilibrium methods and stress transfer range in describing the process of damage. We find that equilibrium approaches are not applicable to the description of damage and the catastrophic failure mechanism if the stress transfer is short ranged. In the long-range limit, equilibrium methods apply only if the healing mechanism associated with ruptured elements is instantaneous. Furthermore we find that the nature of the catastrophic failure depends strongly on the stress transfer range. Long-range transfer systems have a failure mechanism that resembles nucleation. In short-range stress transfer systems, the catastrophic failure is a continuous process that, in some respects, resembles a critical point.

Quantifying uncertainty and sampling quality in biomolecular simulations

Growing computing capacity and algorithmic advances have facilitated the study of increasingly large biomolecular systems at longer timescales. However, with these larger, more complex systems come questions about the quality of sampling and statistical convergence. What size systems can be sampled fully? If a system is not fully sampled, can certain "fast variables" be considered well-converged? How can one determine the statistical significance of observed results? The present review describes statistical tools and the underlying physical ideas necessary to address these questions. Basic definitions and ready-to-use analyses are provided, along with explicit recommendations. Such statistical analyses are of paramount importance in establishing the reliability of simulation data in any given study.

Near-mean-field behavior in the generalized Burridge-Knopoff earthquake model with variable-range stress transfer

Simple models of earthquake faults are important for understanding the mechanisms for their observed behavior in nature, such as Gutenberg-Richter scaling. Because of the importance of long-range interactions in an elastic medium, we generalize the Burridge-Knopoff slider-block model to include variable range stress transfer. We find that the Burridge-Knopoff model with long-range stress transfer exhibits qualitatively different behavior than the corresponding long-range cellular automata models and the usual Burridge-Knopoff model with nearest-neighbor stress transfer, depending on how quickly the friction force weakens with increasing velocity. Extensive simulations of quasiperiodic characteristic events, mode-switching phenomena, ergodicity, and waiting-time distributions are also discussed. Our results are consistent with the existence of a mean-field critical point and have important implications for our understanding of earthquakes and other driven dissipative systems.

Overcoming entropic barrier with coupled sampling at dual resolutions

An enhanced sampling method is proposed for ab initio protein folding simulations. The new method couples a high-resolution model for accuracy and a low-resolution model for efficiency. It aims to overcome the entropic barrier found in the exponentially large protein conformational space when a high-resolution model, such as an all-atom molecular mechanics force field, is used. The proposed method is designed to satisfy the detailed balance condition so that the Boltzmann distribution can be generated in all sampling trajectories in both high and low resolutions. The method was tested on model analytical energy functions and ab initio folding simulations of a beta-hairpin peptide. It was found to be more efficient than replica-exchange method that is used as its building block. Analysis with the analytical energy functions shows that the number of energy calculations required to find global minima and to converge mean potential energies is much fewer with the new method. Ergodic measure shows that the new method explores the conformational space more rapidly. We also studied imperfect low-resolution energy models and found that the introduction of errors in low-resolution models does decrease its sampling efficiency. However, a reasonable increase in efficiency is still observed when the global minima of the low-resolution models are in the vicinity of the global minimum basin of the high-resolution model. Finally, our ab initio folding simulation of the tested peptide shows that the new method is able to fold the peptide in a very short simulation time. The structural distribution generated by the new method at the equilibrium portion of the trajectory resembles that in the equilibrium simulation starting from the crystal structure.

Molecular dynamics study on the solvent dependent heme cooling following ligand photolysis in carbonmonoxy myoglobin

The time scale and mechanism of vibrational energy relaxation of the heme moiety in myoglobin was studied using molecular dynamics simulation. Five different solvent models, including normal water, heavy water, normal glycerol, deuterated glycerol and a nonpolar solvent, and two forms of the heme, one native and one lacking acidic side chains, were studied. Structural alteration of the protein was observed in native myoglobin glycerol solution and native myoglobin water solution. The single-exponential decay of the excess kinetic energy of the heme following ligand photolysis was observed in all systems studied. The relaxation rate depends on the solvent used. However, this dependence cannot be explained using bulk transport properties of the solvent including macroscopic thermal diffusion. The rate and mechanism of heme cooling depends upon the detailed microscopic interaction between the heme and solvent. Three intermolecular energy transfer mechanisms were considered: (i) energy transfer mediated by hydrogen bonds, (ii) direct vibration-vibration energy transfer via resonant interaction, and (iii) energy transfer via vibration-translation or vibration-rotation interaction, or in other words, thermal collision. The hydrogen bond interaction and vibration-vibration interaction between the heme and solvent molecules dominates the energy transfer in native myoglobin aqueous solution and native myoglobin glycerol solutions. For modified myoglobin, the vibration-vibration interaction is also effective in glycerol solution, different from aqueous solution. Thermal collisions form the dominant energy transfer pathway for modified myoglobin in water solution, and for both native myoglobin and modified myoglobin in a nonpolar environment. For native myoglobin in a nonpolar solvent solution, hydrogen bonds between heme isopropionate side chains and nearby protein residues, absent in the modified myoglobin nonpolar solvent solution, are key interactions influencing the relaxation pathways.

Replica exchange with dynamical scaling

A replica exchange method is presented which requires fewer replicas and is designed to be used for large systems. In this method, dynamically scaled replicas are placed between conventional replicas at broadly spaced temperatures. The potential of the scaled replicas is linearly scaled by a dynamical variable which varies between 0 and 1. When the variable is near either end point the replica can undergo exchanges with one of its neighboring replicas. Two different versions of the method are presented for a model system of a small peptide in water. The scaled replica can replace many replicas and the method can be up to ten times more efficient than conventional replica exchange.

Super-Arrhenius diffusion in an undercooled binary Lennard-Jones liquid results from a quantifiable correlation effect

On short time scales an underlying Arrhenius temperature dependence of the diffusion constant can be extracted from the fragile, super-Arrhenius diffusion of a binary Lennard-Jones mixture. This Arrhenius diffusion is related to the true super-Arrhenius behavior by a factor that depends on the average angle between steps in successive time windows. The correction factor accounts for the fact that on average, successive displacements are negatively correlated, and this effect can therefore be linked directly with the higher apparent activation energy for diffusion at low temperature.

Simulation of the Burridge-Knopoff model of earthquakes with variable range stress transfer

Simple models of earthquake faults are important for understanding the mechanisms for their observed behavior, such as Gutenberg-Richter scaling and the relation between large and small events, which is the basis for various forecasting methods. Although cellular automaton models have been studied extensively in the long-range stress transfer limit, this limit has not been studied for the Burridge-Knopoff model, which includes more realistic friction forces and inertia. We find that the latter model with long-range stress transfer exhibits qualitatively different behavior than both the long-range cellular automaton models and the usual Burridge-Knopoff model with nearest-neighbor springs, depending on the nature of the velocity-weakening friction force. These results have important implications for our understanding of earthquakes and other driven dissipative systems.

Diagnosing broken ergodicity using an energy fluctuation metric

The Mountain and Thirumalai energy fluctuation metric, Omega(t), has been used to study the effective ergodicity of 60- and 256-atom binary Lennard-Jones mixtures in order to determine the reliability of the calculated diffusion constants at different energies. A plot of Omega(t) against 1time allows the identification of two distinct regimes: ergodic supercooled liquids, where Omega(t) approaches zero, and nonergodic glassy states, where Omega(t) asymptotically approaches a nonzero value on the molecular dynamics time scale. This approach seems to be more appropriate than attempting to define a threshold value for Omega(t)/Omega(0). The behavior of systems between these two limits, which are nonergodic on the time scale considered but may be approaching ergodicity, was examined for a range of simulation times. The calculated diffusion constants change as effective ergodicity is approached, moving closer to the Vogel-Tammann-Fulcher fit defined by higher-energy systems that are already considered to be effectively ergodic. Using the form of the decay of the metric as a measure of ergodicity, we have been able to reproduce the trend in fragility obtained by Sastry for a 256-atom system [Nature (London) 409, 164 (2001)], correcting some of our earlier results [J. Chem. Phys. 120, 8314 (2004)].

All-exchanges parallel tempering

An alternative exchange strategy for parallel tempering simulations is introduced. Instead of attempting to swap configurations between two randomly chosen but adjacent replicas, the acceptance probabilities of all possible swap moves are calculated a priori. One specific swap move is then selected according to its probability and enforced. The efficiency of the method is illustrated first on the case of two Lennard-Jones (LJ) clusters containing 13 and 31 atoms, respectively. The convergence of the caloric curve is seen to be at least twice as fast as in conventional parallel tempering simulations, especially for the difficult case of LJ31. Further evidence for an improved efficiency is reported on the ergodic measure introduced by Mountain and Thirumalai [J. Phys. Chem. 93, 6975 (1989)], calculated here for LJ13 close to the melting point. Finally, tests on two simple spin systems indicate that the method should be particularly useful when a limited number of replicas are available.

Reentrant behavior of relaxation time with viscosity at varying composition in binary mixtures

In order to understand the long known anomalies in the composition dependence of diffusion and viscosity of binary mixtures, we introduce here two new models and carry out extensive molecular dynamics simulations. In these models, the two molecular species (A and B) have the same diameter and mass. In model I the interspecies interaction is more attractive than that between the pure components, while the reverse is true for model II. Simulations and mode coupling theory calculations reveal that the models can capture a wide variety of behavior observed in experiments, including the reentrant viscosity dependence of relaxation time.

Approach to ergodicity in monte carlo simulations

The approach to the ergodic limit in Monte Carlo simulations is studied using both analytic and numerical methods. With the help of a stochastic model, a metric is defined that enables the examination of a simulation in both the ergodic and nonergodic regimes. In the nonergodic regime, the model implies how the simulation is expected to approach ergodic behavior analytically, and the analytically inferred decay law of the metric allows the monitoring of the onset of ergodic behavior. The metric is related to previously defined measures developed for molecular dynamics simulations, and the metric enables the comparison of the relative efficiencies of different Monte Carlo schemes. Applications to Lennard-Jones 13-particle clusters are shown to match the model for Metropolis, J-walking, and parallel tempering based approaches. The relative efficiencies of these three Monte Carlo approaches are compared, and the decay law is shown to be useful in determining needed high temperature parameters in parallel tempering and J-walking studies of atomic clusters.

Spinodals, scaling, and ergodicity in a threshold model with long-range stress transfer

We present both theoretical and numerical analyses of a cellular automaton version of a slider-block model or threshold model that includes long-range interactions. Theoretically we develop a coarse-grained description in the mean-field (infinite range) limit and discuss the relevance of the metastable state, limit of stability (spinodal), and nucleation to the phenomenology of the model. We also simulate the model and confirm the relevance of the theory for systems with long- but finite-range interactions. Results of particular interest include the existence of Gutenberg-Richter-like scaling consistent with that found on real earthquake fault systems, the association of large events with nucleation near the spinodal, and the result that such systems can be described, in the mean-field limit, with techniques appropriate to systems in equilibrium.

Theoretical probes of conformational fluctuations in S-peptide and RNase A/3'-UMP enzyme product complex

The dynamic properties of the RNase A/3'-UMP enzyme/product complex and the S-peptide of RNase A have been investigated by molecular dynamics simulations using suitable generalization of ideas introduced to probe the energy landscape in structural glasses. We introduce two measures, namely, the kinetic energy fluctuation metric and the force metric, both of which are used to calculate the time needed for sampling the conformation space of the molecules. The calculation of the fluctuation metric requires a single trajectory whereas the force metric is computed using two independent trajectories. The vacuum MD simulations show that for both systems the time required for kinetic energy equipartitioning is surprisingly long even at high temperatures. We show that the force metric is a powerful means of probing the nature and relative importance of conformational substates which determine the dynamics at low temperatures. In particular the time dependence of the non-bonded force metric is used to demonstrate that at low temperatures the system is predominantly localized in a single cluster of conformational substates. The force metric is used to show that relaxation of long range (in sequence space) interactions must be mediated by a sequence of local dihedral angle transitions. We also argue that the time needed for compact structure formation is intimately related to the time needed for the relaxation of the dihedral angle degrees of freedom. The time for non-bonded interactions, which drive protein molecules to fold under appropriate conditions, to relax becomes extremely long as the temperature is lowered suggesting that the formation of maximally compact structure in proteins must be a very slow process.

Exploring the energy landscape in proteins

We present two methods to probe the energy landscape and motions of proteins in the context of molecular dynamics simulations of the helix-forming S-peptide of RNase A and the RNase A-3'-UMP enzyme-product complex. The first method uses the generalized ergodic measure to compute the rate of conformational space sampling. Using the dynamics of nonbonded forces as a means of probing the time scale for ergodicity to be obtained, we argue that even in a relatively short time (< 10 psec) several different conformational substrates are sampled. At longer times, barriers on the order of a few kcal/mol (1 cal = 4.184 J) are involved in the large-scale motion of proteins. We also present an approximate method for evaluating the distribution of barrier heights g(EB) using the instantaneous normal-mode spectra of a protein. For the S-peptide, we show that g(EB) is adequately represented by a Poisson distribution. By comparing with previous work on other systems, we suggest that the statistical characteristics of the energy landscape may be a "universal" feature of all proteins.

The nature of folded states of globular proteins

We suggest, using dynamical simulations of a simple heteropolymer modelling the alpha-carbon sequence in a protein, that generically the folded states of globular proteins correspond to statistically well-defined metastable states. This hypothesis, called the metastability hypothesis, states that there are several free energy minima separated by barriers of various heights such that the folded conformations of a polypeptide chain in each of the minima have similar structural characteristics but have different energies from one another. The calculated structural characteristics, such as bond angle and dihedral angle distribution functions, are assumed to arise from only those configurations belonging to a given minimum. The validity of this hypothesis is illustrated by simulations of a continuum model of a heteropolymer whose low temperature state is a well-defined beta-barrel structure. The simulations were done using a molecular dynamics algorithm (referred to as the "noisy" molecular dynamics method) containing both friction and noise terms. It is shown that for this model there are several distinct metastable minima in which the structural features are similar. Several new methods of analyzing fluctuations in structures belonging to two distinct minima are introduced. The most notable one is a dynamic measure of compactness that can in principle provide the time required for maximal compactness to be achieved. The analysis shows that for a given metastable state in which the protein has a well-defined folded structure the transition to a state of higher compactness occurs very slowly, lending credence to the notion that the system encounters a late barrier in the process of folding to the most compact structure. The examination of the fluctuations in the structures near the unfolding----folding transition temperature indicates that the transition state for the unfolding to folding process occurs closer to the folded state.