Robert K.-Z. Tan

The flying ice cube: Velocity rescaling in molecular dynamics leads to violation of energy equipartition

This article describes an unexpected phenomenon encountered during MD simulations: velocity rescaling using standard protocols can systematically change the proportion of total kinetic energy (KE) found in motions associated with the various degrees of freedom. Under these conditions, the simulation violates the principle of equipartition of energy, which requires a mean kinetic energy of RT/2 in each degree of freedom. A particularly pathological form of this problem occurs if one does not periodically remove the net translation of (and rotation about) the center of mass. In this case, almost all of the kinetic energy is converted into these two kinds of motion, producing a system with almost no kinetic energy associated with the internal degrees of freedom. We call this phenomenon “the flying ice cube.” We present a mathematical analysis of a simple diatomic system with two degrees of freedom, to document the origin of the problem. We then present examples from three kinds of MD simulations, one being an in vacuo simulation on a diatomic system, one involving a low resolution model of DNA in vacuo, and the third using a traditional all‐atom DNA model with full solvation, periodic boundary conditions, and the particle mesh Ewald method for treating long‐range electrostatics. Finally, we discuss methods for avoiding the problem. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 726–740, 1998

Investigation of viral DNA packaging using molecular mechanics models

A simple molecular mechanics model has been used to investigate optimal spool-like packing conformations of double-stranded DNA molecules in viral capsids with icosahedral symmetry. The model represents an elastic segmented chain by using one pseudoatom for each ten basepairs (roughly one turn of the DNA double helix). Force constants for the various terms in the energy function were chosen to approximate known physical properties, and a radial restraint was used to confine the DNA into a sphere with a volume corresponding to that of a typical bacteriophage capsid. When the DNA fills 90% of the spherical volume, optimal packaging is obtained for coaxially spooled models, but this result does not hold when the void volume is larger. When only 60% of the spherical volume is filled with DNA, the lowest energy structure has two layers, with a coiled core packed at an angle to an outer coaxially spooled shell. This relieves bending strain associated with tight curvature near the poles in a model with 100% coaxial spooling. Interestingly, the supercoiling density of these models is very similar to typical values observed in plasmids in bacterial cells. Potential applications of the methodology are also discussed.

Molecular mechanics model of supercoiled DNA

We describe a pseudo-atomic model of supercoiled DNA. Each base-pair of the DNA is represented in the model by three particles placed in a plane. The particle triplets are stacked to model stacked base-pairs in double-helical DNA, and closed circular conformations are generated to investigate supercoiling. This model is less detailed than all-atom models, which are too computationally demanding to be used to study supercoiling. On the other hand, this model contains details at the base-pair level and is therefore more elaborate than elastomechanical models. A potential energy function is written in terms of a set of internal co-ordinates defined to resemble a limited number of helical parameters. The modeled helical parameters, helical twist, base-roll, tilt and rise, are the most important parameters of the global shape of DNA. Experimentally measured mechanical properties of DNA are used to define the forces holding the particles together. We then use a procedure incorporating energy minimization and molecular dynamics to locate low energy conformations of the model DNA. The model was found to behave very much like rubber-tubing and elastomechanical models. The conformations and the effects of supercoiling pressure (a number proportional to the degree to which the total twist of the DNA has been altered from its natural value) on these conformations are all very similar to those observed in the latter two models. We also used this model to examine the effects of supercoiling pressure, base-sequence and mechanical properties on the conformations and energies of five sequences. The sequences studied include models of naturally straight DNA and DNA with static or natural bends.

Modeling large RNAs and ribonucleoprotein particles using molecular mechanics techniques

There is a growing body of low-resolution structural data that can be utilized to devise structural models for large RNAs and ribonucleoproteins. These models are routinely built manually. We introduce an automated refinement protocol to utilize such data for building low-resolution three-dimensional models using the tools of molecular mechanics. In addition to specifying the positions of each nucleotide, the protocol provides quantitative estimates of the uncertainties in those positions, i.e., the resolution of the model. In typical applications, the resolution of the models is about 10-20 A. Our method uses reduced representations and allows us to refine three-dimensional structures of systems as big as the 16S and 23S ribosomal RNAs, which are about one to two orders of magnitude larger than nucleic acids that can be examined by traditional all-atom modeling methods. Nonatomic resolution structural data--secondary structure, chemical cross-links, chemical and enzymatic footprinting patterns, protein positions, solvent accessibility, and so on--are combined with known motifs in RNA structure to predict low-resolution models of large RNAs. These structural constraints are imposed on the RNA chain using molecular mechanics-type potential functions with parameters based on the quality of experimental data. Surface potential functions are used to incorporate shape and positional data from electron microscopy image reconstruction experiments into our models. The structures are optimized using techniques of energy refinement to get RNA folding patterns. In addition to providing a consensus model, the method finds the range of models consistent with the data, which allows quantitative evaluation of the resolution of the model. The method also identifies conflicts in the experimental data. Although our protocol is aimed at much larger RNAs, we illustrate these techniques using the tRNA structure as an example and test-bed.

Yammp: development of a molecular mechanics program using the modular programming method

Molecular mechanics is a fast developing discipline with new methods and potential fields appearing every year. A versatile molecular mechanics program supports many methods and potential fields that make it unavoidably large. There are problems writing and maintaining large programs with traditional methods because of data and other dependencies. Modular programming provides a solution. A program is developed as a collection of highly independent modules containing only related data structures and procedures. These entities are isolated in the module and access to them is provided through a well‐defined and controlled interface. The high degree of independence circumscribes programming errors. Most of all, it reduces the cost of revising the program as only a small part of the program needs to be read and understood for each revision. We implemented a molecular mechanics program, yammp, using the modular programming method.

A Comparison of Six DNA Bending Models

The predictions of six DNA bending models were compared with experimental relative mobility data. The study showed that all the models are reasonably accurate in predicting bending in synthetic sequences and in a natural sequence. The least accurate of these models is the Calladine-Dickerson model. The most consistent model is the ApA Wedge, possibly because it distributes the bends into base-roll and base-tilt components.

Molecular dynamics simulations of small DNA plasmids: Effects of sequence and supercoiling on intramolecular motions

Small (600 base pair) DNA plasmids were modeled with a simplified representation (3DNA) and the intramolecular motions were studied using molecular mechanics and molecular dynamics techniques. The model is detailed enough to incorporate sequence effects. At the same time, it is simple enough to allow long molecular dynamics simulations. The simulations revealed that large-scale slithering occurs in a homogeneous sequence. In a heterogeneous sequence, containing numerous small intrinsic curves, the centers of the curves are preferentially positioned at the tips of loops. With more curves than loop tips (two in unbranched supercoiled DNA), the heterogeneous sequence plasmid slithers short distances to reposition other curves into the loop tips. However, the DNA is immobilized most of the time, with the loop tips positioned over a few favored curve centers. Branching or looping also appears in the heterogeneous sequence as a new method of repositioning the loop tips. Instead of a smooth progression of increasing writhing with increasing linking difference, theoretical studies have predicted that there is a threshold between unwrithed and writhed DNA at a linking difference between one and two. This has previously been observed in simulations of static structures and is demonstrated here for dynamic homogeneous closed DNA. Such an abrupt transition is not found in the heterogeneous sequence in both the static and dynamic cases.

MOLECULAR MODELING OF CLOSED CIRCULAR DNA THERMODYNAMIC ENSEMBLES

Many modeling studies of supercoiled DNA are based on equilibrium structures from theoretical calculations or energy minimization. Since closed circular DNAs are flexible, it is possible that errors are introduced by calculating properties from a single minimum energy structure, rather than from a complete thermodynamic ensemble. We have investigated this question using molecular dynamics simulations on a low resolution molecular mechanics model in which each base pair is represented by three points (a plane). This allows the inclusion of sequence-dependent variations of tip, inclination, and twist. Three kinds of sequences were tested: (1) homogeneous DNA, in which all base pairs have the helicoidal parameters of an ideal, average B-DNA; (2) random sequence DNA; and (3) curved DNA. We examined the rate of convergence of various structural parameters. Convergence for most of these is slowest for homogeneous sequences, more rapid for random sequences, and most rapid for curved sequences. The most slowly converging parameter is the antipodes profile. In a plasmid with N base pairs (bp), the antipodes distance is the distance dij from base pair i to base pair j halfway around the plasmid, j = i + N/2. The antipodes profile at time tau is a plot of dij over the range i = 1, N/2. In a homogeneous plasmid, convergence requires that the antipodes profile averaged over time must be flat. Even in the small plasmids examined here, the average properties of the ensembles were found to differ from those of static equilibrium structures. These effects will be even more dramatic for larger plasmids. Further, average and dynamic properties are affected by both plasmid size and sequence.

DNA topological context affects access to eukaryotic DNA topoisomerase I.

We have analyzed the reactivity of a 217 base pair segment of the intrinsically curved Crithidia fasciculata kinetoplast DNA towards eukaryotic DNA topoisomerase I. The substrates were open [linear fragment and nicked circle] and closed minidomains [closed relaxed circle and circles with linking differences of -1 and -2]. We interpreted the results with the aid of a model that was used to predict the structures of the topoisomers. The modelling shows that the delta Lk(-1) form is unusually compact because of the curvature in the DNA. To determine the role of sequence-directed curvature in both the experimental and modeling studies, controls were examined in which the curved Crithidia sequence was replaced by an uncurved sequence obtained from the plasmid pBR322. Reactivity of the Crithidia DNA [as analyzed both by the cleavage and topoisomerization reactions] markedly varied among the DNA forms: (i) the hierarchy of overall reactivity observed is: linear fragment > nicked circular, closed circular [delta Lk(0)], interwound [delta Lk(-2)] > bent interwound [delta Lk(-1)]; (ii) the intensity of several cleavage positions differs among DNA forms. The results show that eukaryotic DNA topoisomerase I is very sensitive to the conformation of the substrates and that its reactivity is modulated by the variation of the compactness of the DNA molecule. The C. fasciculata sequence contains a highly curved segment that determines the conformation of the closed circle in a complex way.

AUGUR: a program to predict, display and analyze the tertiary structure of B-DNA

AUGUR is a program to predict, display and analyze the three-dimensional structure of B-DNA. The user can choose one of six models to predict the helical parameters of a given sequence. These parameters are then used to generate the coordinates of the DNA model in three-dimensional space (trajectory). The trajectory can be displayed and rotated on a graphics terminal. The trajectory and helical parameters can also be searched for bends and structural homologues.