Extracting parameters for base-pair level models of DNA from molecular dynamics simulations

 A method is described to extract a complete set of sequence-dependent energy parameters for a rigid base-pair model of DNA from molecular dynamics (MD) simulations. The method is properly consistent with equilibrium statistical mechanics and leads to effective inertia parameters for the base-pair units as well as stacking and stiffness parameters for the base-pair junctions. We give explicit formulas that yield a complete set of base-pair model parameters in terms of equilibrium averages that can be estimated from a time series generated in an MD simulation. The expressions to be averaged depend strongly both on the choice of coordinates used to describe rigid-body orientations and on the choice of strain measures at each junction. Received: 12 July 2000 / Accepted: 5 January 2001 / Published online: 3 May 2001     

ConCept: de novo design of synthetic receptors for targeted ligands

Low-molecular-weight receptors that bind targeted guest molecules have a wide range of potential applications but are difficult to design. This paper describes an evolutionary method for computer-aided design of such receptors that works by linking together chemical components from a user-defined library around a stable conformation of the targeted ligand. The software can operate in three modes: de novo design, in which it builds a wide variety of receptors from small components; macrocycle design, in which it builds homopolymeric macrocycles around the ligand; and elaboration of an existing receptor structure. The top candidates generated by the automatic construction process are further studied with detailed affinity calculations whose validity is supported by prior studies of experimentally characterized host-guest systems. All three modes of operation are illustrated here through the design of novel adenine receptors.     

Extracting atoms from molecular electron densities via integral equations

The observation that a molecular electron density is close to the superposition of its constituent atoms leads naturally to the idea of modeling a density by a sum of nuclear-centered, spherically symmetric functions. The functions that are optimal in a least-squares sense are known as Stewart atoms. Previous attempts to construct Stewart atoms by expanding them in an auxiliary basis have been thwarted by slow convergence with respect to the size of the auxiliary basis used. We present a method for constructing Stewart atoms via convolution integrals which bypasses the need for an auxiliary basis, and is able to produce highly accurate approximations to Stewart atoms.     

Cause and effect reversed in non-equilibrium molecular dynamics: an easy route to transport coefficients

A novel non-equilibrium method for calculating transport coefficients is presented. It reverses the experimental cause-and-effect picture, e.g. for the calculation of viscosities: the effect, the momentum flux or stress, is imposed, whereas the cause, the velocity gradient or shear rates, is obtained from the simulation. It differs from other Norton-ensemble methods by the way, in which the steady-state fluxes are maintained. This method involves a simple exchange of particle momenta, which is easy to implement and to analyse. Moreover, it can be made to conserve the total energy as well as the total linear momentum, so no thermostatting is needed. The resulting raw data are robust and rapidly converging. The method is tested on the calculation of the shear viscosity, the thermal conductivity and the Soret coefficient (thermal diffusion) for the Lennard–Jones (LJ) fluid near its triple point. Possible applications to other transport coefficients and more complicated systems are discussed.     

Polar solvation dynamics of lysozyme from molecular dynamics studies

The solvation dynamics of a protein are believed to be sensitive to its secondary structures. We have explored such sensitivity in this article by performing room temperature molecular dynamics simulation of an aqueous solution of lysozyme. Nonuniform long-time relaxation patterns of the solvation time correlation function for different segments of the protein have been observed. It is found that relatively slower long-time solvation components of the α-helices and β-sheets of the protein are correlated with lower exposure of their polar probe residues to bulk solvent and hence stronger interactions with the dynamically restricted surface water molecules. These findings can be verified by appropriate experimental studies.     

Estimating entropies from molecular dynamics simulations

While the determination of free-energy differences by MD simulation has become a standard procedure for which many techniques have been developed, total entropies and entropy differences are still hardly ever computed. An overview of techniques to determine entropy differences is given, and the accuracy and convergence behavior of five methods based on thermodynamic integration and perturbation techniques was evaluated using liquid water as a test system. Reasonably accurate entropy differences are obtained through thermodynamic integration in which many copies of a solute are desolvated. When only one solute molecule is involved, only two methods seem to yield useful results, the calculation of solute-solvent entropy through thermodynamic integration, and the calculation of solvation entropy through the temperature derivative of the corresponding free-energy difference. One-step perturbation methods seem unsuitable to obtain entropy estimates.     

Vesicle shapes from molecular dynamics simulations

Lipid bilayer membranes are known to form various structures such as large sheets or vesicles. When the two leaflets of the bilayer have an equal composition, the membrane preferentially forms a flat sheet or a spherical vesicle. However, a difference in the composition of the two leaflets may result in a curved bilayer or in a wide variety of vesicle shapes. Vesicles with different shapes have already been shown in experiments and diverse vesicle shapes have been predicted theoretically from energy minimization of continuous curves. Here we present a molecular dynamics study of the effect of small changes in the phospholipid headgroups on the spontaneous curvature of the bilayer and on the resulting vesicle shape transformations. Small asymmetries in the bilayers already result in high spontaneous curvature and large vesicle deformations. Vesicle shapes that are formed include ellipsoids, discoids, pear-shaped vesicles, cup-shaped vesicles, as well as budded vesicles. Comparison of these vesicles with theoretically derived vesicle shapes shows both resemblances and differences.     

Extracting effective normal modes from equilibrium dynamics at finite temperature

A general method for obtaining effective normal modes of a molecular system from molecular dynamics simulations is presented. The method is based on a localization criterion for the Fourier transformed velocity time-correlation functions of the effective modes. For a given choice of the localization function used, the method becomes equivalent to the principal mode analysis (PMA) based on covariance matrix diagonalization. On the other hand, a proper choice of the localization function leads to a novel method with a strong analogy with the usual normal mode analysis of equilibrium structures, where the Hessian system at the minimum energy structure is replaced by the thermal averaged Hessian, although the Hessian itself is never actually calculated. This method does not introduce any extra numerical cost during the simulation and bears the same simplicity as PMA itself. It can thus be readily applied to ab initio molecular dynamics simulations. Three such examples are provided here. First we recover effective normal modes of an isolated formaldehyde molecule computed at 20 K in very good agreement with the results of a normal mode analysis performed at its equilibrium structure. We then illustrate the applicability of the method for liquid phase studies. The effective normal modes of a water molecule in liquid water and of a uracil molecule in aqueous solution can be extracted from ab initio molecular dynamics simulations of these two systems at 300 K.     

Solution dynamics and molecular structure elucidation of novel aluminium derivatives containing diaminodimethylsilane ligands

The interaction of dimethyldiaminosilane ligands of general formula SiMe2(NR2)(NR'2)(NR2, NR'2 = NiHPr, NHtBu, NC4H8, NHCH2CH2NMe2) with AlX3 (X = Cl, Me) has been investigated and the synthesis of novel aluminium derivatives is reported, namely AlMe3[SiMe2(NR2)(NR'2)], AlX2[SiMe2(NR)(NR'2)] and AlMe[SiMe2(NR)2], containing the silane ligand as neutral, monoanionic and dianionic species, respectively. Moreover, the solution molecular structures and dynamics have been elucidated via 1D/2D variable temperature NMR spectroscopy showing the influence of the N-substituents of the silane ligand and of the aluminium ancillary ligands.     

Extracting elements of molecular structure from the all-particle wave function

Structural information is extracted from the all-particle (non-Born-Oppenheimer) wave function by calculating radial and angular densities derived from n-particle densities. As a result, one- and two-dimensional motifs of classical molecular structure can be recognized in quantum mechanics. Numerical examples are presented for three- (H(-), Ps(-), H(2)(+)), four- (Ps(2), H(2)), and five-particle (H(2)D(+)) systems.     

Characterization of polyethylene crystallization from an oriented melt by molecular dynamics simulation

Molecular dynamics is used to characterize the process of crystallization for a united atom model of polyethylene. An oriented melt is produced by uniaxial deformation under constant load, followed by quenching below the melting temperature at zero load. The development of crystallinity is monitored simultaneously using molecular-based order parameters for density, energy, and orientation. For crystallization temperatures ranging from 325 to 375 K, these simulations clearly show the hallmarks of crystal nucleation and growth. We can identify multiple nucleation events, lamellar growth up to the limit imposed by periodic boundaries of the simulation cell, and lamellar thickening. We observe a competition between the rate of nucleation, which results in multiple crystallites, the rate of chain extension, which results in thicker lamellae, and the rate of chain conformational relaxation, which is manifested in lower degrees of residual order in the noncrystalline portion of the simulation. The temperature dependence of lamellar thickness is in accord with experimental data. At the higher temperatures, tilted chain lamellae are observed to form with lamellar interfaces corresponding approximately to the [201] facet, indicative of the influence of interfacial energy.     

Measuring coexisting densities from a two-phase molecular dynamics simulation by voronoi tessellations

A new algorithm is presented that allows for the determination of bulk liquid and vapor densities from a two-phase Molecular Dynamics (2phiMD) simulation. This new method does not use any arbitrary cutoffs for phase definitions; rather it uses single-phase simulations as a self-consistency check. The method does not use any spatial bins for generating histograms of local properties, thereby avoiding the statistical issues associated with bins. Finally, it allows one to approach very close to the critical point. The new method utilizes Voronoi tessellations to determine the molecular volume of every point at every instance in a molecular dynamics simulation. Since the molecular volume is calculated throughout the simulation, statistical parameters such as the average molecular volume and average molecular variance are easy to obtain. To define the phases, the normalized variance of the molecular volume from 1phiMD and 2phiMD is used as a self-consistency check. The new method gives new insight into the nature of the near-subcritical fluid. The critical properties from this analysis are T(c) = 1.293 and rho(c) = 0.313. Direct simulation of the two-phase system was performed up to a temperature of 1.292. The results show excellent agreement to experimental results and Gibbs Ensemble Monte Carlo for coexisting densities. We see that well below the critical temperature, some particles are neither liquid nor vapor. These interfacial particles are primarily, but not exclusively, concentrated at the bulk interface. However, as we approach the critical point, some particles are considered both liquid and vapor. These interfacial particles are distributed through the system.     

Multivalent cyclooligosaccharides: versatile carbohydrate clusters with dual role as molecular receptors and lectin ligands

Results obtained over the past decade towards the preparation of multitopic carbohydrate architectures combining the molecular inclusion capabilities of cyclomaltooligosaccharide receptors (cyclodextrins, CDs) and the recognition properties of saccharide ligands towards biological receptors are discussed. The potential of these new sugar-based "intelligent" transporters for site specific delivery of therapeutics is outlined.     

Thermal conductivities of molecular liquids by reverse nonequilibrium molecular dynamics

The reverse nonequilibrium molecular dynamics method for thermal conductivities is adapted to the investigation of molecular fluids. The method generates a heat flux through the system by suitably exchanging velocities of particles located in different regions. From the resulting temperature gradient, the thermal conductivity is then calculated. Different variants of the algorithm and their combinations with other system parameters are tested: exchange of atomic velocities versus exchange of molecular center-of-mass velocities, different exchange frequencies, molecular models with bond constraints versus models with flexible bonds, united-atom versus all-atom models, and presence versus absence of a thermostat. To help establish the range of applicability, the algorithm is tested on different models of benzene, cyclohexane, water, and n-hexane. We find that the algorithm is robust and that the calculated thermal conductivities are insensitive to variations in its control parameters. The force field, in contrast, has a major influence on the value of the thermal conductivity. While calculated and experimental thermal conductivities fall into the same order of magnitude, in most cases the calculated values are systematically larger. United-atom force fields seem to do better than all-atom force fields, possibly because they remove high-frequency degrees of freedom from the simulation, which, in nature, are quantum-mechanical oscillators in their ground state and do not contribute to heat conduction.     

Extracting the time scales of conformational dynamics from single-molecule single-photon fluorescence statistics

The quenching rate of a fluorophore attached to a macromolecule can be rather sensitive to its conformational state. The decay of the corresponding fluorescence lifetime autocorrelation function can therefore provide unique information on the time scales of conformational dynamics. The conventional way of measuring the fluorescence lifetime autocorrelation function involves evaluating it from the distribution of delay times between photoexcitation and photon emission. However, the time resolution of this procedure is limited by the time window required for collecting enough photons in order to establish this distribution with sufficient signal-to-noise ratio. Yang and Xie have recently proposed an approach for improving the time resolution, which is based on the argument that the autocorrelation function of the delay time between photoexcitation and photon emission is proportional to the autocorrelation function of the square of the fluorescence lifetime [Yang, H.; Xie, X. S. J. Chem. Phys. 2002, 117, 10965]. In this paper, we show that the delay-time autocorrelation function is equal to the autocorrelation function of the square of the fluorescence lifetime divided by the autocorrelation function of the fluorescence lifetime. We examine the conditions under which the delay-time autocorrelation function is approximately proportional to the autocorrelation function of the square of the fluorescence lifetime. We also investigate the correlation between the decay of the delay-time autocorrelation function and the time scales of conformational dynamics. The results are demonstrated via applications to a two-state model and an off-lattice model of a polypeptide.     

Micellar crystals in solution from molecular dynamics simulations

Polymers with both soluble and insoluble blocks typically self-assemble into micelles, which are aggregates of a finite number of polymers where the soluble blocks shield the insoluble ones from contact with the solvent. Upon increasing concentration, these micelles often form gels that exhibit crystalline order in many systems. In this paper, we present a study of both the dynamics and the equilibrium properties of micellar crystals of triblock polymers using molecular dynamics simulations. Our results show that equilibration of single micelle degrees of freedom and crystal formation occur by polymer transfer between micelles, a process that is described by transition state theory. Near the disordered (or melting) transition, bcc lattices are favored for all triblocks studied. Lattices with fcc ordering are also found but only at lower kinetic temperatures and for triblocks with short hydrophilic blocks. Our results lead to a number of theoretical considerations and suggest a range of implications to experimental systems with a particular emphasis on Pluronic polymers.     

Optimized explicit-solvent replica exchange molecular dynamics from scratch

Replica exchange molecular dynamics (REMD) simulations have become an important tool to study proteins and other biological molecules in silico. However, such investigations require considerable, and often prohibitive, numerical effort when the molecules are simulated in explicit solvents. In this communication we show that in this case the cost can be minimized by choosing the number of replicas as N(opt) approximately 1+0.594 radical C ln(Tmax/Tmin), where C is the specific heat, and the temperatures distributed according to Ti(opt) approximately T min(Tmax/Tmin)(i-1)/(N-1).