Sensitivity of molecular dynamics simulations to the choice of the X-ray structure used to model an enzymatic reaction

A subject of great practical importance that has not received much attention is the question of the sensitivity of molecular dynamics simulations to the initial X-ray structure used to set up the calculation. We have found two cases in which seemingly similar structures lead to quite different results, and in this article we present a detailed analysis of these cases. The first case is acyl-CoA dehydrogenase, and the chief difference of the two structures is attributed to a slight shift in a backbone carbonyl that causes a key residue (the proton-abstracting base) to be in a bad conformation for reaction. The second case is xylose isomerase, and the chief difference of the two structures appears to be the ligand sphere of a Mg2+ metal cofactor that plays an active role in catalysis.     

Genetic assignment methods for the direct, real-time estimation of migration rate: a simulation-based exploration of accuracy and power

Genetic assignment methods use genotype likelihoods to draw inference about where individuals were or were not born, potentially allowing direct, real-time estimates of dispersal. We used simulated data sets to test the power and accuracy of Monte Carlo resampling methods in generating statistical thresholds for identifying F0 immigrants in populations with ongoing gene flow, and hence for providing direct, real-time estimates of migration rates. The identification of accurate critical values required that resampling methods preserved the linkage disequilibrium deriving from recent generations of immigrants and reflected the sampling variance present in the data set being analysed. A novel Monte Carlo resampling method taking into account these aspects was proposed and its efficiency was evaluated. Power and error were relatively insensitive to the frequency assumed for missing alleles. Power to identify F0 immigrants was improved by using large sample size (up to about 50 individuals) and by sampling all populations from which migrants may have originated. A combination of plotting genotype likelihoods and calculating mean genotype likelihood ratios (DLR) appeared to be an effective way to predict whether F0 immigrants could be identified for a particular pair of populations using a given set of markers.     

Sequence alignments and pair hidden Markov models using evolutionary history

This work presents a novel pairwise statistical alignment method based on an explicit evolutionary model of insertions and deletions (indels). Indel events of any length are possible according to a geometric distribution. The geometric distribution parameter, the indel rate, and the evolutionary time are all maximum likelihood estimated from the sequences being aligned. Probability calculations are done using a pair hidden Markov model (HMM) with transition probabilities calculated from the indel parameters. Equations for the transition probabilities make the pair HMM closely approximate the specified indel model. The method provides an optimal alignment, its likelihood, the likelihood of all possible alignments, and the reliability of individual alignment regions. Human alpha and beta-hemoglobin sequences are aligned, as an illustration of the potential utility of this pair HMM approach.     

Simulating structural and thermodynamic properties of carcinogen-damaged DNA

A pair of stereoisomeric covalent adducts to guanine in double-stranded DNA, derived from the reaction of mutagenic and tumorigenic metabolites of benzo[a]pyrene, have been well characterized structurally and thermodynamically. Both high-resolution NMR solution structures and an array of thermodynamic data are available for these 10S (+)- and 10R (-)-trans-anti -[BP]-N(2)-dG adducts in double-stranded deoxyoligonucleotides. The availability of experimentally well-characterized duplexes containing these two stereoisomeric guanine adducts provides an opportunity for evaluating the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method for computing thermodynamic properties from molecular dynamics ensembles. We have carried out 3-ns molecular dynamics simulations, using NMR solution structures as the starting models for the 10S (+)- and 10R (-)-trans-anti-dG adducts in a DNA duplex 11-mer using AMBER 6.0. We employed the MM-PBSA method to compute the free energies, enthalpies, and entropies of the two adducts. Our complete thermodynamic analysis agrees quite well with the full experimental thermodynamic characterization of these adducts, showing essentially equal stabilities of the two adducts. We also calculated the nuclear Overhauser effect (NOE) distances from the molecular dynamics trajectories, and compared them against the experimental NMR-derived NOE distances. Our results showed that the simulated structures are in good agreement with the NMR experimental NOE data. Furthermore, the molecular dynamics simulations provided new structural and biological insights. Specifically, the puzzling observation that the BP aromatic ring system in the 10S (+)-trans-anti-dG adduct is more exposed to the aqueous solvent than the 10R (-)-trans-anti-dG adduct, is rationalized in terms of the adduct structures. The structural and thermodynamic features of these stereoisomeric adducts are also discussed in relation to their reported low susceptibilities to nucleotide excision repair.     

Sequence context dependence of tandem guanine:adenine mismatch conformations in RNA: a continuum solvent analysis

Guanine:adenine (G:A) mismatches and in particular tandem G:A (tG:A) mismatches are frequently observed in biological RNA molecules and can serve as sites for tertiary interaction, metal binding and protein recognition. Depending on the surrounding sequence tG:A mismatches can adopt different basepairing topologies. In the sequence context (5'-) GGAC (tandem G:A in bold) a face-to-face (imino or Watson-Crick-like) pairing is preferred whereas in the CGAG context, G and A adopt a sheared arrangement. Systematic conformational searches with a generalized Born continuum model and molecular dynamics simulations including explicit water molecules and ions have been used to generate face-to-face and sheared tG:A mismatches in both CGAG and GGAC sequence contexts. Conformations from both approaches were evaluated using the same force field and a Poisson-Boltzmann continuum solvent model. Although the substate analysis predicted the sheared arrangement to be energetically preferred in both sequence contexts, a significantly greater preference of the sheared form was found for the CGAG context. In agreement with the experimental observation, the analysis of molecular dynamics trajectories indicated a preference of the sheared form in the case of the CGAG-context and a favorization of the face-to-face form in the case of the GGAC context. The computational studies allowed to identify energetic contributions that stabilize or destabilize the face-to-face and sheared tandem mismatch topologies. The calculated nonpolar solvation and Lennard-Jones packing interaction were found to stabilize the sheared topology independent of the sequence context. Electrostatic contributions are predicted to make the most significant contribution to the sequence context dependence on the structural preference of tG:A mismatches.     

Conformational behavior of chondroitin and chondroitin sulfate in relation to their physical properties as inferred by molecular modeling

Chondroitin and chondroitin sulfates belong to the family of glycosaminoglycans. They are most widely distributed in animal tissues, where they are involved in structural functions and in cell-cell communication. Their basic structures consist of a disaccharidic repeating unit of beta-D-glucuronic acid (GlcA) and 2-acetamido-2-deoxy-beta-D-galactose (GalNAc), this latter being sulfated at different positions. Molecular mechanics has been applied to calculate the adiabatic energy maps for each of the constituting disaccharides of chondroitin, chondroitin 4-sulfate, and chondroitin 6-sulfate using the MM3 force field. Based on these maps, higher levels of structural organization have been simulated. On one hand, the disordered state is studied through a Metropolis-based algorithm; the resulting chains present a behavior of semirigid polymers, with an order of stiffness: chondroitin 4-sulfate > chondroitin > chondroitin 6-sulfate. On the other hand, the exploration of the stable ordered forms leads to numerous helical conformations of comparable energies. Several of these conformations correspond to the experimentally observed ones. The ability of coordination with cations has also been explored, resulting in a preferential stereospecificity for calcium ions when compared to sodium ions.     

Predicting the conformational states of cyclic tetrapeptides

Biologically active cyclic tetrapeptides, usually found among fungi metabolites, exhibit phytotoxic or cytostatic activities that are likely to be governed by specific conformations adopted in solution. For conformational studies and drug design, there is a strong interest in using fast and reliable methods to determine correctly the conformational population of cyclotetrapeptides. We show here that standard molecular mechanics computational approach gives satisfactory results. The method was validated step by step by experimental data either obtained after synthesis and NMR analysis, or found in the literature. The cyclo(Gly)(4), cyclo(Ala)(4), cyclo(Sar)(4), and cyclo(SarGly)(2) peptides were used to evaluate the prediction of the peptide backbone conformation, and the detailed conformational analysis of tentoxin, a natural phytotoxic cyclotetrapeptide in which N-alkylated peptide bonds alternate with regular secondary ones, was used to validate the computation of conformers proportions. From the knowledge of an initial cyclic primary structure and of the D or L configuration of the amino acids, we show that it is possible to determine the exact orientation of carbonyl groups and to predict the nature of conformers present in solution. The proportion of each conformer can be inferred from a statistical thermodynamics approach by using the potential energy values of each conformer, computed by molecular mechanics methods with the TRIPOS force field, which allowed us to account for the solvent. The solvent contribution was processed by two different methods according to the nature of the interactions: whether through the dielectric constant introduced in the electrostatic potential, when interaction with solute molecules are weak or negligible, or through the computation of free energy of solvation using the algorithm SILVERWARE for solvents explicitly interacting with the solute. When applied to tentoxin, this conformational analysis yielded results in very good agreement with the experimental data reported by Pinet et al. (Biopolymers, 1995, Vol. 36, pp. 135-152), on both the nature of existing conformers and their relative proportions, whatever the nature of the considered solvent.     

Molecular dynamics simulations of an oligosaccharide using a force field modified for carbohydrates

Parameterization of the phi and omega torsion angles in pyranosidic saccharides was performed based on density functional theory calculations. The modified CHARMM force field, which is referred to as PARM22/SU01, was tested on a glucosyl trisaccharide. A molecular dynamics simulation of the oligosaccharide with explicit water as solvent was performed to investigate the conformational flexibility. Protonz.sbnd;proton distances and heteronuclear spin-spin coupling constants were calculated from the trajectories and showed good agreement to those previously determined by NMR spectroscopy.     

Toward a realistic model of mutations affecting fitness

Analysis of a recent mutation accumulation (MA) experiment has led to the suggestion that as many as one-half of spontaneous mutations in Arabidopsis are advantageous for fitness. We evaluate this in the light of data from other MA experiments, along with molecular evidence, that suggest the vast majority of new mutations are deleterious.