Recurrence quantification analysis as a tool for the characterization of molecular dynamics simulations

A molecular dynamics simulation of a Lennard-Jones fluid, and a trajectory of the B1 immunoglobulin G-binding domain of streptococcal protein G (B1-IgG) simulated in water are analyzed by recurrence quantification, which is noteworthy for its independence from stationarity constraints, as well as its ability to detect transients, and both linear and nonlinear state changes. The results demonstrate the sensitivity of the technique for the discrimination of phase sensitive dynamics. Physical interpretation of the recurrence measures is also discussed.Comment: 7 pages, 8 figures, revtex; revised for review for Phys. Rev. E (clarifications and expansion of discussion)-- addition of the 8 postscript figures previously omitted, but unchanged from version

Molecular dynamics simulation of polymer helix formation using rigid-link methods

Molecular dynamics simulations are used to study structure formation in simple model polymer chains that are subject to excluded volume and torsional interactions. The changing conformations exhibited by chains of different lengths under gradual cooling are followed until each reaches a state from which no further change is possible. The interactions are chosen so that the true ground state is a helix, and a high proportion of simulation runs succeed in reaching this state; the fraction that manage to form defect-free helices is a function of both chain length and cooling rate. In order to demonstrate behavior analogous to the formation of protein tertiary structure, additional attractive interactions are introduced into the model, leading to the appearance of aligned, antiparallel helix pairs. The simulations employ a computational approach that deals directly with the internal coordinates in a recursive manner; this representation is able to maintain constant bond lengths and angles without the necessity of treating them as an algebraic constraint problem supplementary to the equations of motion.Comment: 15 pages, 14 figure

Maximum Flux Transition Paths of Conformational Change

Given two metastable states A and B of a biomolecular system, the problem is to calculate the likely paths of the transition from A to B. Such a calculation is more informative and more manageable if done for a reduced set of collective variables chosen so that paths cluster in collective variable space. The computational task becomes that of computing the "center" of such a cluster. A good way to define the center employs the concept of a committor, whose value at a point in collective variable space is the probability that a trajectory at that point will reach B before A. The committor "foliates" the transition region into a set of isocommittors. The maximum flux transition path is defined as a path that crosses each isocommittor at a point which (locally) has the highest crossing rate of distinct reactive trajectories. (This path is different from that of the MaxFlux method of Huo and Straub.) It is argued that such a path is nearer to an ideal path than others that have been proposed with the possible exception of the finite-temperature string method path. To make the calculation tractable, three approximations are introduced, yielding a path that is the solution of a nonsingular two-point boundary-value problem. For such a problem, one can construct a simple and robust algorithm. One such algorithm and its performance is discussed.Comment: 7 figure

Integrity of H1 helix in prion protein revealed by molecular dynamic simulations to be especially vulnerable to changes in the relative orientation of H1 and its S1 flank

In the template-assistance model, normal prion protein (PrPC), the pathogenic cause of prion diseases such as Creutzfeldt-Jakob (CJD) in human, Bovine Spongiform Encephalopathy (BSE) in cow, and scrapie in sheep, converts to infectious prion (PrPSc) through an autocatalytic process triggered by a transient interaction between PrPC and PrPSc. Conventional studies suggest the S1-H1-S2 region in PrPC to be the template of S1-S2 $\beta$-sheet in PrPSc, and the conformational conversion of PrPC into PrPSc may involve an unfolding of H1 in PrPC and its refolding into the $\beta$-sheet in PrPSc. Here we conduct a series of simulation experiments to test the idea of transient interaction of the template-assistance model. We find that the integrity of H1 in PrPC is vulnerable to a transient interaction that alters the native dihedral angles at residue Asn$^{143}$, which connects the S1 flank to H1, but not to interactions that alter the internal structure of the S1 flank, nor to those that alter the relative orientation between H1 and the S2 flank.Comment: A major revision on statistical analysis method has been made. The paper now has 23 pages, 11 figures. This work was presented at 2006 APS March meeting session K29.0004 at Baltimore, MD, USA 3/13-17, 2006. This paper has been accepted for pubcliation in European Biophysical Journal on Feb 2, 200

Isomorphs in model molecular liquids

Isomorphs are curves in the phase diagram along which a number of static and dynamic quantities are invariant in reduced units. A liquid has good isomorphs if and only if it is strongly correlating, i.e., the equilibrium virial/potential energy fluctuations are more than 90% correlated in the NVT ensemble. This paper generalizes isomorphs to liquids composed of rigid molecules and study the isomorphs of two systems of small rigid molecules, the asymmetric dumbbell model and the Lewis-Wahnstrom OTP model. In particular, for both systems we find that the isochoric heat capacity, the excess entropy, the reduced molecular center-of-mass self part of the intermediate scattering function, the reduced molecular center-of-mass radial distribution function to a good approximation are invariant along an isomorph. In agreement with theory, we also find that an instantaneous change of temperature and density from an equilibrated state point to another isomorphic state point leads to no relaxation. The isomorphs of the Lewis-Wahnstrom OTP model were found to be more approximative than those of the asymmetric dumbbell model, which is consistent with the OTP model being less strongly correlating. For both models we find "master isomorphs", i.e., isomorphs have identical shape in the virial/potential energy phase diagram.Comment: 20 page

Growth of (110) Diamond using pure Dicarbon

We use a density-functional based tight-binding method to study diamond growth steps by depositing dicarbon species onto a hydrogen-free diamond (110) surface. Subsequent C_2 molecules are deposited on an initially clean surface, in the vicinity of a growing adsorbate cluster, and finally, near vacancies just before completion of a full new monolayer. The preferred growth stages arise from C_2n clusters in near ideal lattice positions forming zigzag chains running along the [-110] direction parallel to the surface. The adsorption energies are consistently exothermic by 8--10 eV per C_2, depending on the size of the cluster. The deposition barriers for these processes are in the range of 0.0--0.6 eV. For deposition sites above C_2n clusters the adsorption energies are smaller by 3 eV, but diffusion to more stable positions is feasible. We also perform simulations of the diffusion of C_2 molecules on the surface in the vicinity of existing adsorbate clusters using an augmented Lagrangian penalty method. We find migration barriers in excess of 3 eV on the clean surface, and 0.6--1.0 eV on top of graphene-like adsorbates. The barrier heights and pathways indicate that the growth from gaseous dicarbons proceeds either by direct adsorption onto clean sites or after migration on top of the existing C_2n chains.Comment: 8 Pages, 7 figure

Coarse-Graining Protein Structures With Local Multivariate Features from Molecular Dynamics

A multivariate statistical theory, local feature analysis (LFA), extracts functionally relevant domains from molecular dynamics (MD) trajectories. The LFA representations, like those of principal component analysis (PCA), are low dimensional and provide a reduced basis set for collective motions of simulated proteins, but the local features are sparsely distributed and spatially localized, in contrast to global PCA modes. One key problem in the assignment of local features is the coarse-graining of redundant LFA output functions by means of seed atoms. One can solve the combinatorial problem by adding seed atoms one after another to a growing set, minimizing a reconstruction error at each addition. This allows for an efficient implementation, but the sequential algorithm does not guarantee the optimal mutual correlation of the sequentially assigned features. Here, we present a novel coarse-graining algorithm for proteins that directly minimizes the mutual correlation of seed atoms by Monte Carlo (MC) simulations. Tests on MD trajectories of two biological systems, bacteriophage T4 lysozyme and myosin II motor domain S1, demonstrate that the new algorithm provides statistically reproducible results and describes functionally relevant dynamics. The well-known undersampling of large-scale motion by short MD simulations is apparent also in our model, but the new coarse-graining offers a major advantage over PCA; converged features are invariant across multiple windows of the trajectory, dividing the protein into converged regions and a smaller number of localized, undersampled regions. In addition to its use in structure classification, the proposed coarse-graining thus provides a localized measure of MD sampling efficiency

A critical assessment of methods for the intrinsic analysis of liquid interfaces. 1. surface site distributions

Substantial progress in our understanding of interfacial structure and dynamics has stemmed from the recent development of algorithms that allow for an intrinsic analysis of fluid interfaces. These work by identifying the instantaneous location of the interface, at the atomic level, for each molecular configuration and then computing properties relative to this location. Such a procedure eliminates the broadening of the interface caused by capillary waves and reveals the underlying features of the system. However, a precise definition of which molecules actually belong to the interfacial layer is difficult to achieve in practice. Furthermore, it is not known if the different intrinsic analysis methods are consistent with each other and yield similar results for the interfacial properties. In this paper, we carry out a systematic and detailed comparison of the available methods for intrinsic analysis of fluid interfaces, based on a molecular dynamics simulation of the interface between liquid water and carbon tetrachloride. We critically assess the advantages and shortcomings of each method, based on reliability, robustness, and speed of computation, and establish consistent criteria for determining which molecules belong to the surface layer. We believe this will significantly contribute to make intrinsic analysis methods widely and routinely applicable to interfacial systems

Protecting High Energy Barriers: A New Equation to Regulate Boost Energy in Accelerated Molecular Dynamics Simulations

Molecular dynamics (MD) is one of the most common tools in computational chemistry. Recently, our group has employed accelerated molecular dynamics (aMD) to improve the conformational sampling over conventional molecular dynamics techniques. In the original aMD implementation, sampling is greatly improved by raising energy wells below a predefined energy level. Recently, our group presented an alternative aMD implementation where simulations are accelerated by lowering energy barriers of the potential energy surface. When coupled with thermodynamic integration simulations, this implementation showed very promising results. However, when applied to large systems, such as proteins, the simulation tends to be biased to high energy regions of the potential landscape. The reason for this behavior lies in the boost equation used since the highest energy barriers are dramatically more affected than the lower ones. To address this issue, in this work, we present a new boost equation that prevents oversampling of unfavorable high energy conformational states. The new boost potential provides not only better recovery of statistics throughout the simulation but also enhanced sampling of statistically relevant regions in explicit solvent MD simulations

Differential Base Stacking Interactions Induced by Trimethylene Interstrand DNA Cross-Links in the 5′-CpG-3′ and 5′-GpC-3′ Sequence Contexts

Synthetically derived trimethylene interstrand DNA cross-links have been used as surrogates for the native cross-links that arise from the 1,N 2-deoxyguanosine adducts derived from R,β-unsaturated aldehydes. The native enal-mediated cross-linking occurs in the 5′-CpG-3 ′ sequence context but not in the 5′-GpC-3 ′ sequence context. The ability of the native enal-derived 1,N 2-dG adducts to induce interstrand DNA cross-links in the 5′-CpG-3 ′ sequence as opposed to the 5′-GpC-3 ′ sequence is attributed to the destabilization of the DNA duplex in the latter sequence context. Here, we report higher accuracy solution structures of the synthetically derived trimethylene cross-links, which are refined from NMR data with the AMBER force field. When the synthetic trimethylene cross-links are placed into either the 5′-CpG-3′ or the 5′-GpC-3 ′ sequence contexts, the DNA duplex maintains B-DNA geometry with structural perturbations confined to the cross-linked base pairs. Watson-Crick hydrogen bonding is conserved throughout the duplexes. Although different from canonical B-DNA stacking, the cross-linked and the neighbor base pairs stack in the 5′-CpG-3 ′ sequence. In contrast, the stacking at the cross-linked base pairs in the 5′-GpC-3 ′ sequence is greatly perturbed. The π-stacking interactions between the crosslinked and the neighbor base pairs are reduced. This is consistent with remarkable chemical shift perturbations of the C 5 H5 and H6 nucleobase protons that shifted downfield by 0.4-0.5 ppm. In contrast