Molecular dynamics simulations of the 136 unique tetranucleotide sequences of DNA oligonucleotides. I. Research design and results on d(CpG) steps

MD-DATA: the legacy of the ABC Consortium

The ABC Consortium has been generating nucleic-acids MD trajectories for more than 20 years. This brief comment highlights the importance of this data for the field, which triggered a number of critical studies, including force-field parameterization and development of new coarse-grained and mesoscopic models. With the world entering into a new data-driven era led by artificial intelligence, where data is becoming more essential than ever, the ABC initiative is leading the way for nucleic acid flexibility.

Multiscale simulations of DNA from electrons to nucleosomes: 22 years of the Ascona B-DNA Consortium

This commentary provides a retrospective on the Ascona B-DNA Consortium (ABC) initiative and on the conference held in April 2023 at Ascona, Switzerland, where we celebrated 22 years of the consortium, sharing the latest advances in simulations and experiments of the effects of sequence on the mechanical properties of DNA from electrons to nucleosomes.

Systematic assessment of the flexibility of uracil damaged DNA

Uracil is a common DNA lesion which is recognized and removed by uracil DNA-glycosylase (UDG) as a part of the base excision repair pathway. Excision proceeds by base flipping, and UDG efficiency is thought to depend on the ease of deformability of the bases neighboring the lesion. We used molecular dynamics simulations to assess the flexibility of a large library of dsDNA strands, containing all tetranucleotide motifs with U:A, U:G, T:A or C:G base pairs. Our study demonstrates that uracil damaged DNA largely follows trends in flexibility of undamaged DNA. Measured bending persistence lengths, groove widths, step parameters and base flipping propensities demonstrate that uracil increases the flexibility of DNA, and that U:G base paired strands are more flexible than U:A strands. Certain sequence contexts are more deformable than others, with a key role for the 3' base next to uracil. Flexibilities are large when this base is an A or G, and repressed for a C or T. A 5' T adjacent to the uracil strongly promotes flexibility, but other 5' bases are less influential. DNA bending is correlated to step deformations and base flipping, and bending aids flipping. Our study implies that the link between substrate flexibility and UDG efficiency is widely valid, helps explain why UDG prefers to bind U:G base paired strands, and suggests that the DNA bending angle of the UDG-substrate complex is optimal for base flipping.Communicated by Ramaswamy H. Sarma.

From complex data to clear insights: visualizing molecular dynamics trajectories

Advances in simulations, combined with technological developments in high-performance computing, have made it possible to produce a physically accurate dynamic representation of complex biological systems involving millions to billions of atoms over increasingly long simulation times. The analysis of these computed simulations is crucial, involving the interpretation of structural and dynamic data to gain insights into the underlying biological processes. However, this analysis becomes increasingly challenging due to the complexity of the generated systems with a large number of individual runs, ranging from hundreds to thousands of trajectories. This massive increase in raw simulation data creates additional processing and visualization challenges. Effective visualization techniques play a vital role in facilitating the analysis and interpretation of molecular dynamics simulations. In this paper, we focus mainly on the techniques and tools that can be used for visualization of molecular dynamics simulations, among which we highlight the few approaches used specifically for this purpose, discussing their advantages and limitations, and addressing the future challenges of molecular dynamics visualization.

Sequence-Dependent Properties of the RNA Duplex

Sequence-dependent properties of the DNA duplex have been accurately described using extensive molecular dynamics simulations. The RNA duplex meanwhile─which is typically represented as a sequence-averaged rigid rod─does not benefit from having equivalent molecular dynamics simulations. In this paper, we present a massive simulation effort using a set of ABC-optimized duplexes from which we derived tetramer-resolution properties of the RNA duplex and a simple mesoscopic model that can represent elastic properties of long RNA duplexes. Despite the extreme chemical similarity between DNA and RNA, the local and global elastic properties of the duplexes are very different. DNA duplexes show a complex and nonelastic pattern of flexibility, for instance, while RNA duplexes behave as an elastic system whose deformations can be represented by simple harmonic potentials. In RNA duplexes (RNA2), not only are intra- and interbase pair parameters (equilibrium and mechanical) different from those in the equivalent DNA duplex sequences (DNA2) but the correlations between movements also differ. Simple statements on the relative flexibility or stability of both polymers are meaningless and should be substituted by a more detailed description depending on the sequence and the type of deformation considered.

Molecular dynamics simulation-based trinucleotide and tetranucleotide level structural and energy characterization of the functional units of genomic DNA

Genomes of most organisms on earth are written in a universal language of life, made up of four units - adenine (A), thymine (T), guanine (G), and cytosine (C), and understanding the way they are put together has been a great challenge to date. Multiple efforts have been made to annotate this wonderfully engineered string of DNA using different methods but they lack a universal character. In this article, we have investigated the structural and energetic profiles of both prokaryotes and eukaryotes by considering two essential genomic sites, viz., the transcription start sites (TSS) and exon-intron boundaries. We have characterized these sites by mapping the structural and energy features of DNA obtained from molecular dynamics simulations, which considers all possible trinucleotide and tetranucleotide steps. For DNA, these physicochemical properties show distinct signatures at the TSS and intron-exon boundaries. Our results firmly convey the idea that DNA uses the same dialect for prokaryotes and eukaryotes and that it is worth going beyond sequence-level analyses to physicochemical space to determine the functional destiny of DNA sequences.

Revisiting DNA Sequence-Dependent Deformability in High-Resolution Structures: Effects of Flanking Base Pairs on Dinucleotide Morphology and Global Chain Configuration

DNA carries more than the list of biochemical ingredients that drive the basic functions of living systems. The sequence of base pairs includes a multitude of structural and energetic signals, which determine the degree to which the long, threadlike molecule moves and how it responds to proteins and other molecules that control its processing and govern its packaging. The chemical composition of base pairs directs the spatial disposition and fluctuations of successive residues. The observed arrangements of these moieties in high-resolution protein–DNA crystal structures provide one of the best available estimates of the natural, sequence-dependent structure and deformability of the double-helical molecule. Here, we update the set of knowledge-based elastic potentials designed to describe the observed equilibrium structures and configurational fluctuations of the ten unique base-pair steps. The large number of currently available structures makes it possible to characterize the configurational preferences of the DNA base-pair steps within the context of their immediate neighbors, i.e., tetrameric context. Use of these knowledge-based potentials shows promise in accounting for known effects of sequence in long chain molecules, e.g., the degree of curvature reported in classic gel mobility studies and the recently reported sequence-dependent responses of supercoiled minicircles to nuclease cleavage.

emDNA - A Tool for Modeling Protein-decorated DNA Loops and Minicircles at the Base-pair Step Level

Computational modeling of nucleic acids plays an important role in molecular biology, enhancing our general understanding of the relationship between structure and function. Biophysical studies have provided a wealth of information on how double-helical DNA responds to proteins and other molecules in its local environment but far less understanding of the larger scale structural responses found in protein-decorated loops and minicircles. Current computational models of DNA range from detailed all-atom molecular dynamics studies, which produce rich time and spatially dependent depictions of small DNA fragments, to coarse-grained simulations, which sacrifice detailed physical and chemical information to treat larger scale systems. The treatment of DNA used here, at the base-pair step level with rigid-body parameters, allows one to develop quality models hundreds of base pairs long from local, sequence-specific features found from experiment. The emDNA software takes advantage of this framework, producing elastically optimized models of DNA at thermal equilibrium with built-in or user-generated elastic models. This versatile program, in combination with case studies included in this article, allows users of any skill level to develop and investigate mesoscale models of their own design. The functionality of emDNA includes a tool to incorporate experiment-specific configurations, e.g. protein-bound and/or melted DNA from known high-resolution structures, within higher-order 3D models by fixing the orientation and position of user-specified base pairs. The software provides a new avenue into multiscale genetic modeling, giving a wide range of users a deeper understanding of DNA mesoscale organization and the opportunity to pose new questions in genetic research. The publicly available emDNA software, including build instructions and usage information, is available on GitHub (https://nicocvn.github.io/emDNA/).

Computational investigation of the impact of core sequence on immobile DNA four-way junction structure and dynamics

Immobile four-way junctions (4WJs) are core structural motifs employed in the design of programmed DNA assemblies. Understanding the impact of sequence on their equilibrium structure and flexibility is important to informing the design of complex DNA architectures. While core junction sequence is known to impact the preferences for the two possible isomeric states that junctions reside in, previous investigations have not quantified these preferences based on molecular-level interactions. Here, we use all-atom molecular dynamics simulations to investigate base-pair level structure and dynamics of four-way junctions, using the canonical Seeman J1 junction as a reference. Comparison of J1 with equivalent single-crossover topologies and isolated nicked duplexes reveal conformational impact of the double-crossover motif. We additionally contrast J1 with a second junction core sequence termed J24, with equal thermodynamic preference for each isomeric configuration. Analyses of the base-pair degrees of freedom for each system, free energy calculations, and reduced-coordinate sampling of the 4WJ isomers reveal the significant impact base sequence has on local structure, isomer bias, and global junction dynamics.

Systematic Comparison of the Structural and Dynamic Properties of Commonly Used Water Models for Molecular Dynamics Simulations

Water is a unique solvent that is ubiquitous in biology and present in a variety of solutions, mixtures, and materials settings. It therefore forms the basis for all molecular dynamics simulations of biological phenomena, as well as for many chemical, industrial, and materials investigations. Over the years, many water models have been developed, and it remains a challenge to find a single water model that accurately reproduces all experimental properties of water simultaneously. Here, we report a comprehensive comparison of structural and dynamic properties of 30 commonly used 3-point, 4-point, 5-point, and polarizable water models simulated using consistent settings and analysis methods. For the properties of density, coordination number, surface tension, dielectric constant, self-diffusion coefficient, and solvation free energy of methane, models published within the past two decades consistently show better agreement with experimental values compared to models published earlier, albeit with some notable exceptions. However, no single model reproduced all experimental values exactly, highlighting the need to carefully choose a water model for a particular study, depending on the phenomena of interest. Finally, machine learning algorithms quantified the relationship between the water model force field parameters and the resulting bulk properties, providing insight into the parameter-property relationship and illustrating the challenges of developing a water model that can accurately reproduce all properties of water simultaneously.

Effects of Structural Dynamics on Charge Carrier Transfer in B-DNA: A Combined MD and RT-TDDFT Study

Hole transfer along the axis of duplex DNA has been the focus of physical chemistry research for decades, with implications in diverse fields, from nanotechnology to cell oxidative damage. Computational approaches are particularly amenable for this problem, to complement experimental data for interpretation of transfer mechanisms. To be predictive, computational results need to account for the inherent mobility of biological molecules during the time frame of experimental measurements. Here, we address the structural variability of B-DNA and its effects on hole transfer in a combined molecular dynamics (MD) and real-time time-dependent density functional theory (RT-TDDFT) study. Our results show that quantities that characterize the charge transfer process, such as the time-dependent dipole moment and hole population at a specific site, are sensitive to structural changes that occur on the nanosecond time scale. We extend the range of physical properties for which such a correlation has been observed, further establishing the fact that quantitative computational data on charge transfer properties should include statistical averages. Furthermore, we use the RT-TDDFT results to assess an efficient tight-binding method suitable for high-throughput predictions. We demonstrate that charge transfer, although affected by structural variability, on average, remains strong in AA and GG dimers.

Intron exon boundary junctions in human genome have in-built unique structural and energetic signals

Precise identification of correct exon–intron boundaries is a prerequisite to analyze the location and structure of genes. The existing framework for genomic signals, delineating exon and introns in a genomic segment, seems insufficient, predominantly due to poor sequence consensus as well as limitations of training on available experimental data sets. We present here a novel concept for characterizing exon–intron boundaries in genomic segments on the basis of structural and energetic properties. We analyzed boundary junctions on both sides of all the exons (3 28 368) of protein coding genes from human genome (GENCODE database) using 28 structural and three energy parameters. Study of sequence conservation at these sites shows very poor consensus. It is observed that DNA adopts a unique structural and energy state at the boundary junctions. Also, signals are somewhat different for housekeeping and tissue specific genes. Clustering of 31 parameters into four derived vectors gives some additional insights into the physical mechanisms involved in this biological process. Sites of structural and energy signals correlate well to the positions playing important roles in pre-mRNA splicing.

A novel method SEProm for prokaryotic promoter prediction based on DNA structure and energetics

MOTIVATION

Despite conservation in general architecture of promoters and protein-DNA interaction interface of RNA polymerases among various prokaryotes, identification of promoter regions in the whole genome sequences remains a daunting challenge. The available tools for promoter prediction do not seem to address the problem satisfactorily, apparently because the biochemical nature of promoter signals is yet to be understood fully. Using 28 structural and 3 energetic parameters, we found that prokaryotic promoter regions have a unique structural and energy state, quite distinct from that of coding regions and the information for this signature state is in-built in their sequences. We developed a novel promoter prediction tool from these 31 parameters using various statistical techniques.

RESULTS

Here, we introduce SEProm, a novel tool that is developed by studying and utilizing the in-built structural and energy information of DNA sequences, which is applicable to all prokaryotes including archaea. Compared to five most recent, diverged and current best available tools, SEProm performs much better, predicting promoters with an 'F-value' of 82.04 and 'Precision' of 81.08. The next best 'F-value' was obtained with PromPredict (72.14) followed by BProm (68.37). On the basis of 'Precision' value, the next best 'Precision' was observed for Pepper (75.39) followed by PromPredict (72.01). SEProm maintained the lead even when comparison was done on two test organisms (not involved in training for SEProm).

AVAILABILITY AND IMPLEMENTATION

The software is freely available with easy to follow instructions (www.scfbio-iitd.res.in/software/TSS_Predict.jsp).

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Atomistic insight into sequence-directed DNA bending and minicircle formation propensity in the absence and presence of phased A-tracts

Bending of double-stranded (ds) DNA plays a crucial role in many important biological processes and is relevant for nanotechnological applications. Among all the elements that have been studied in relation to dsDNA bending, A-tracts stand out as one of the most controversial. The "ApA wedge" theory was disproved when a series of linear polynucleotides containing phased 5'-A₄T₄-3' or 5'-T₄A₄-3' runs were shown to be bent or straight, respectively, and crystallographic evidence revealed that A-tracts are unbent. Furthermore, some of the smallest dsDNA minicircles described to date (~ 100 bp in size) lack A-tracts and are subjected to varying levels of torsional stress. Representative DNA sequences from this experimental background were modeled in atomic detail and their dynamic behavior was simulated over hundreds of nanoseconds using the AMBER force field ParmBSC1. Subsequent analysis of the resulting trajectories allowed us to (i) unambiguously establish the location of the bends in all cases; (ii) identify the structural elements that are directly responsible for the macroscopically detected curvature; and (iii) reveal the importance not only of coherently summing the effects of the bending elements when they are in synchrony with the natural repeat of the helix (i.e. separated by an integral number of helical turns) but also when alternated with a half-integral separation of opposite effects. We conclude that the major determinant of the macroscopically observed bending is the proper grouping and phasing of the positive roll imposed by pyrimidine-purine (YR) steps and the negative or null roll characteristic of RY steps and A-tracts, respectively. This conclusion is in very good agreement with the structural parameters experimentally derived for much smaller DNA molecules either alone or as found in DNA-protein complexes. We expect that this work will pave the way for future studies on drug-induced DNA bending, DNA shape readout by transcription factors, structure of circular extrachromosomal DNA, and custom design of curved DNA origami scaffolds.

The static and dynamic structural heterogeneities of B-DNA: extending Calladine-Dickerson rules

We present a multi-laboratory effort to describe the structural and dynamical properties of duplex B-DNA under physiological conditions. By processing a large amount of atomistic molecular dynamics simulations, we determine the sequence-dependent structural properties of DNA as expressed in the equilibrium distribution of its stochastic dynamics. Our analysis includes a study of first and second moments of the equilibrium distribution, which can be accurately captured by a harmonic model, but with nonlocal sequence-dependence. We characterize the sequence-dependent choreography of backbone and base movements modulating the non-Gaussian or anharmonic effects manifested in the higher moments of the dynamics of the duplex when sampling the equilibrium distribution. Contrary to prior assumptions, such anharmonic deformations are not rare in DNA and can play a significant role in determining DNA conformation within complexes. Polymorphisms in helical geometries are particularly prevalent for certain tetranucleotide sequence contexts and are always coupled to a complex network of coordinated changes in the backbone. The analysis of our simulations, which contain instances of all tetranucleotide sequences, allow us to extend Calladine-Dickerson rules used for decades to interpret the average geometry of DNA, leading to a set of rules with quantitative predictive power that encompass nonlocal sequence-dependence and anharmonic fluctuations.

Modulation of the helical properties of DNA: next-to-nearest neighbour effects and beyond

We used extensive molecular dynamics simulations to study the structural and dynamic properties of the central d(TpA) step in the highly polymorphic d(CpTpApG) tetranucleotide. Contrary to the assumption of the dinucleotide-model and its nearest neighbours (tetranucleotide-model), the properties of the central d(TpA) step change quite significantly dependent on the next-to-nearest (hexanucleotide) sequence context and in a few cases are modulated by even remote neighbours (beyond next-to-nearest from the central TpA). Our results highlight the existence of previously undescribed dynamical mechanisms for the transmission of structural information into the DNA and demonstrate the existence of certain sequences with special physical properties that can impact on the global DNA structure and dynamics.

TMB Library of Nucleosome Simulations

Nucleosomes are the fundamental building blocks of chromatin, the biomaterial that houses the genome in all higher organisms. A nucleosome consists of 145-147 base pairs of DNA wrapped 1.7 times around eight histones. Given a four-letter code (A, C, G, T), there are approximately 4147 or 1088 oligonucleotides that can form a nucleosome. Comparative, rather than comprehensive, studies are required. Here we introduce the TMB Library of nucleosome simulations and present a meta-analysis of over 20 μs of all atom molecular dynamics simulations representing 518 different realizations of the nucleosome. The TMB Library serves as a reference for future comparative, on-demand simulations of nucleosomes and a demonstration of iBIOMES Lite as a tool for managing a laboratory's simulation library. For every simulation, dewatered trajectories, RMSD, and DNA helical parameter data are provided through iBIOMES Lite in a Web browser and a file browser format. A novel view of nucleosomal DNA emerges from our meta-analysis of the TMB Library. DNA conformation is restricted to a specific left-handed superhelix, but the range of conformations observed for individual bases and base pairs is not more restricted nor more highly deformed than DNA free in solution. With the exception of Roll, mean DNA helical parameter values obtained from simulations of nucleosomes are largely within the range of thermal motion of DNA free in solution. The library provides evidence of DNA kinking in the nucleosome and clearly demonstrates the effects of DNA sequence on the gross structure and dynamics of nucleosomes. These effects and mispositioning of the 601 super strong nucleosome positioning sequence can be detected in short simulations (10 ns). Collectively, the results provide a basis for comparative simulation studies of nucleosomes and extend our understanding of the binding of proteins and drugs to nucleosomal DNA. The TMB Library can be found at http://dna.engr.latech.edu/~tmbshare/ .

Mechanistic role of nucleobases in self-cleavage catalysis of hairpin ribozyme at ambient versus high-pressure conditions

Ribozymes catalyze the site-specific self-cleavage of intramolecular phosphodiester bonds. Initially thought to act as metalloenzymes, they are now known to be functional even in the absence of divalent metal ions and specific nucleobases directly participate in the self-cleavage reaction. Here, we use extensive replica exchange molecular dynamics simulations to probe the precise mechanistic role of nucleobases by simulating precatalytic reactant and active precursor states of a hairpin ribozyme along its reaction path at ambient as well as high-pressure conditions. The results provide novel key insights into the self-cleavage of ribozymes. We find that deprotonation of the hydroxyl group is crucial and might be the penultimate step to the self-cleavage. The G8 nucleobase is found to stabilize the activated precursor into inline arrangement for facile nucleophilic attack of the scissile phosphate only after deprotonation of the hydroxyl group. The protonated A38 nucleobase, in contrast, mainly acts a proton donor to the O5'-oxygen leaving group that eventually leads to the self-cleavage. Indeed, systematic high-pressure simulations of catalytically relevant states confirm these findings and, moreover, provide support to the role of ribozymes as piezophilic biocatalysts with regard to their relevance in early life under extreme conditions in the realm of RNA world hypothesis.

C5' omitted DNA enhances bendability and protein binding

Protein-DNA interactions are of great biological importance. The specificity and strength of these intimate contacts are crucial in the proper functioning of a cell, wherein the role of DNA dynamic bendability has been a matter of discussion. We relate DNA bendability to protein binding by introducing some simple modifications in the DNA structure. We removed C5' carbon in first modified structure and the second has an additional carbon between C3' and 3'-OH, hereby pronounced as C(-) and C(+) nucleic acids respectively. We observed that C(+) nucleic acid retains B-DNA duplex as seen by means of 500 ns long molecular dynamics (MD) simulations, structural and energetic calculations, while C(-) nucleic acid attains a highly bend structure. We transferred these observations to a protein-DNA system in order to monitor as to what extent the bendability enhances the protein binding. The energetics of binding is explored by performing 100 ns long MD simulations on control and modified DNA-protein complexes followed by running MM-PBSA/GBSA calculations on the resultant structures. It is observed that C(+) nucleic acid has protein binding in close correspondence to the control system (∼-14 kcal/mol) due to their relatable structure, while the C(-) nucleic acid displayed high binding to the protein (∼-18 kcal/mol). DelPhi based calculations reveal that the high binding could be the result of enhanced electrostatic interactions caused by exposed bases in the bend structure for protein recognition. Such modified oligonucleotides, due to their improved binding to protein and resistance to nuclease degradation, have a great therapeutic value.

ChemGenome2.1: An Ab Initio Gene Prediction Software

Gene prediction, also known as gene identification, gene finding, gene recognition, or gene discovery, is among one of the important problems of molecular biology and is receiving increasing attention due to the advent of large-scale genome sequencing projects. We designed an ab initio model (called ChemGenome) for gene prediction in prokaryotic genomes based on physicochemical characteristics of codons. In this chapter, we present the methodology of the latest version of this model ChemGenome2.1 (CG2.1). The first module of the protocol builds a three-dimensional vector from three calculated quantities for each codon-the double-helical trinucleotide base pairing energy, the base pair stacking energy, and an index of the propensity of a codon for protein-nucleic acid interactions. As this three-dimensional vector moves along any genome, the net orientation of the resultant vector should differ significantly for gene and non-genic regions to make a distinction feasible. The predicted putative protein-coding genes from above parameters are passed through a second module of the protocol which reduces the number of false positives by utilizing a filter based on stereochemical properties of protein sequences. The chemical properties of amino acid side chains taken into consideration are the presence of sp3 hybridized γ carbon atom, hydrogen bond donor ability, short/absence of δ carbon and linearity of the side chains/non-occurrence of bi-dentate forks with terminal hydrogen atoms in the side chain. The final prediction of the potential protein-coding genes is based on the frequency of occurrence of amino acids in the predicted protein sequences and their deviation from the frequency values of Swissprot protein sequences, both at monomer and tripeptide levels. The final screening is based on Z-score. Though CG2.1 is a gene finding tool for prokaryotes, considering the underlying similarity in the chemical and physical properties of DNA among prokaryotes and eukaryotes, we attempted to evaluate its applicability for gene finding in the lower eukaryotes. The results give a hope that the concept of gene finding based on physicochemical model of codons is a viable idea for eukaryotes as well, though, undoubtedly, improvements are needed.

Physicochemical property based computational scheme for classifying DNA sequence elements of Saccharomyces cerevisiae

GenerationE of huge "omics" data necessitates the development and application of computational methods to annotate the data in terms of biological features. In the context of DNA sequence, it is important to unravel the hidden physicochemical signatures. For this purpose, we have considered various sequence elements such as promoter, ACS, LTRs, telomere, and retrotransposon of the model organism Saccharomyces cerevisiae. Contributions due to di-nucleotides play a major role in studying the DNA conformation profile. The physicochemical parameters used are hydrogen bonding energy, stacking energy and solvation energy per base pair. Our computational study shows that all sequence elements in this study have distinctive physicochemical signatures and the same can be exploited for prediction experiments. The order that we see in a DNA sequence is dictated by biological regions and hence, there exists role of dependency in the sequence makeup, keeping this in mind we are proposing two computational schemes (a) using a windowing block size procedure and (b) using di-nucleotide transitions. We obtained better discriminating profile when we analyzed the sequence data in windowing manner. In the second novel approach, we introduced the di-nucleotide transition probability matrix (DTPM) to study the hidden layer of information embedded in the sequences. DTPM has been used as weights for scanning and predictions. This proposed computational scheme incorporates the memory property which is more realistic to study the physicochemical properties embedded in DNA sequences. Our analysis shows that the DTPM scheme performs better than the existing method in this applied region. Characterization of these elements will be a key to genome editing applications and advanced machine learning approaches may also require such distinctive profiles as useful input features.

Toward a Universal Structural and Energetic Model for Prokaryotic Promoters

With almost no consensus promoter sequence in prokaryotes, recruitment of RNA polymerase (RNAP) to precise transcriptional start sites (TSSs) has remained an unsolved puzzle. Uncovering the underlying mechanism is critical for understanding the principle of gene regulation. We attempted to search the hidden code in ∼16,500 promoters of 12 prokaryotes representing two kingdoms in their structure and energetics. Twenty-eight fundamental parameters of DNA structure including backbone angles, basepair axis, and interbasepair and intrabasepair parameters were used, and information was extracted from x-ray crystallography data. Three parameters (solvation energy, hydrogen-bond energy, and stacking energy) were selected for creating energetics profiles using in-house programs. DNA of promoter regions was found to be inherently designed to undergo a change in every parameter undertaken for the study, in all prokaryotes. The change starts from some distance upstream of TSSs and continues past some distance from TSS, hence giving a signature state to promoter regions. These signature states might be the universal hidden codes recognized by RNAP. This observation was reiterated when randomly selected promoter sequences (with little sequence conservation) were subjected to structure generation; all developed into very similar three-dimensional structures quite distinct from those of conventional B-DNA and coding sequences. Fine structural details at important motifs (viz. -11, -35, and -75 positions relative to TSS) of promoters reveal novel to our knowledge and pointed insights for RNAP interaction at these locations; it could be correlated with how some particular structural changes at the -11 region may allow insertion of RNAP amino acids in interbasepair space as well as facilitate the flipping out of bases from the DNA duplex.

Sequence-Dependent Three Interaction Site Model for Single- and Double-Stranded DNA

We develop a robust coarse-grained model for single- and double-stranded DNA by representing each nucleotide by three interaction sites (TIS) located at the centers of mass of sugar, phosphate, and base. The resulting TIS model includes base-stacking, hydrogen bond, and electrostatic interactions as well as bond-stretching and bond angle potentials that account for the polymeric nature of DNA. The choices of force constants for stretching and the bending potentials were guided by a Boltzmann inversion procedure using a large representative set of DNA structures extracted from the Protein Data Bank. Some of the parameters in the stacking interactions were calculated using a learning procedure, which ensured that the experimentally measured melting temperatures of dimers are faithfully reproduced. Without any further adjustments, the calculations based on the TIS model reproduce the experimentally measured salt and sequence-dependence of the size of single-stranded DNA (ssDNA), as well as the persistence lengths of poly(dA) and poly(dT) chains. Interestingly, upon application of mechanical force, the extension of poly(dA) exhibits a plateau, which we trace to the formation of stacked helical domains. In contrast, the force-extension curve (FEC) of poly(dT) is entropic in origin and could be described by a standard polymer model. We also show that the persistence length of double-stranded DNA, formed from two complementary ssDNAs, is consistent with the prediction based on the worm-like chain. The persistence length, which decreases with increasing salt concentration, is in accord with the Odijk-Skolnick-Fixman theory intended for stiff polyelectrolyte chains near the rod limit. Our model predicts the melting temperatures of DNA hairpins with excellent accuracy, and we are able to recover the experimentally known sequence-specific trends. The range of applications, which did not require adjusting any parameter after the initial construction based solely on PDB structures and melting profiles of dimers, attests to the transferability and robustness of the TIS model for ssDNA and dsDNA.

Physics-based all-atom modeling of RNA energetics and structure

The database of RNA sequences is exploding, but knowledge of energetics, structures, and dynamics lags behind. All-atom computational methods, such as molecular dynamics, hold promise for closing this gap. New algorithms and faster computers have accelerated progress in improving the reliability and accuracy of predictions. Currently, the methods can facilitate refinement of experimentally determined nuclear magnetic resonance and x-ray structures, but are 'unreliable' for predictions based only on sequence. Much remains to be discovered, however, about the many molecular interactions driving RNA folding and the best way to approximate them quantitatively. The large number of parameters required means that a wide variety of experimental results will be required to benchmark force fields and different approaches. As computational methods become more reliable and accessible, they will be used by an increasing number of biologists, much as x-ray crystallography has expanded. Thus, many fundamental physical principles underlying the computational methods are described. This review presents a summary of the current state of molecular dynamics as applied to RNA. It is designed to be helpful to students, postdoctoral fellows, and faculty who are considering or starting computational studies of RNA. WIREs RNA 2017, 8:e1422. doi: 10.1002/wrna.1422.

Theoretical characterization of the conformational features of unnatural oligonucleotides containing a six nucleotide genetic alphabet

The addition of the unnatural P:Z base pair to the four naturally occurring DNA bases expands the genetic alphabet and yields an artificially expanded genetic information system (AEGIS). Herein, the structural feature of oligonucleotides containing a novel unnatural P:Z base pair is characterized using both molecular dynamics and quantum chemistry. The results show that the incorporation of the novel artificial base pair (P:Z) preserves the global conformational feature of duplex DNA except for some local structures. The Z-nitro group imparts new properties to the groove width, which widens the major groove. The unnatural oligonucleotides containing mismatched base pairs exhibit low stability. This ensures efficient and high-fidelity replication. In general, the incorporation of the P:Z pair strengthens the stability of the corresponding DNA duplex. The calculated results also show that the thermostability originates from both hydrogen interaction and stacking interaction. The Z-nitro group plays an important role in enhancing the stability of the H-bonds and stacking strength of the P:Z pair. Overall, the present results provide theoretical insights in the exploration of artificially expanded genetic information systems.

Theoretical insights into the structural and fluorescence properties of DNA containing fluorescent nucleobases

The structural and fluorescence properties of DNA containing fluorescent nucleobases are explored in detail. In particular, the effect of stacking with different neighboring bases on fluorescence properties is elucidated.

Physico-chemical fingerprinting of RNA genes

We advance here a novel concept for characterizing different classes of RNA genes on the basis of physico-chemical properties of DNA sequences. As knowledge-based approaches could yield unsatisfactory outcomes due to limitations of training on available experimental data sets, alternative approaches that utilize properties intrinsic to DNA are needed to supplement training based methods and to eventually provide molecular insights into genome organization. Based on a comprehensive series of molecular dynamics simulations of Ascona B-DNA consortium, we extracted hydrogen bonding, stacking and solvation energies of all combinations of DNA sequences at the dinucleotide level and calculated these properties for different types of RNA genes. Considering ∼7.3 million mRNA, 255 524 tRNA, 40 649 rRNA (different subunits) and 5250 miRNA, 3747 snRNA, gene sequences from 9282 complete genome chromosomes of all prokaryotes and eukaryotes available at NCBI, we observed that physico-chemical properties of different functional units on genomic DNA differ in their signatures.

How accurate are accurate force-fields for B-DNA?

Last generation of force-fields are raising expectations on the quality of molecular dynamics (MD) simulations of DNA, as well as to the belief that theoretical models can substitute experimental ones in several cases. However these claims are based on limited benchmarks, where MD simulations have shown the ability to reproduce already existing ‘experimental models’, which in turn, have an unclear accuracy to represent DNA conformation in solution. In this work we explore the ability of different force-fields to predict the structure of two new B-DNA dodecamers, determined herein by means of ¹H nuclear magnetic resonance (NMR). The study allowed us to check directly for experimental NMR observables on duplexes previously not solved, and also to assess the reliability of ‘experimental structures’. We observed that technical details in the annealing procedures can induce non-negligible local changes in the final structures. We also found that while not all theoretical simulations are equally reliable, those obtained using last generation of AMBER force-fields (BSC1 and BSC0_OL15) show predictive power in the multi-microsecond timescale and can be safely used to reproduce global structure of DNA duplexes and fine sequence-dependent details.

Analyzing DNA curvature and its impact on the ionic environment: application to molecular dynamics simulations of minicircles

We propose a method for analyzing the magnitude and direction of curvature within nucleic acids, based on the curvilinear helical axis calculated by Curves+. The method is applied to analyzing curvature within minicircles constructed with varying degrees of over- or under-twisting. Using the molecular dynamics trajectories of three different minicircles, we are able to quantify how curvature varies locally both in space and in time. We also analyze how curvature influences the local environment of the minicircles, notably via increased heterogeneity in the ionic distributions surrounding the double helix. The approach we propose has been integrated into Curves+ and the utilities Canal (time trajectory analysis) and Canion (environmental analysis) and can be used to study a wide variety of static and dynamic structural data on nucleic acids.

[KCl] Dependence of B-DNA Groove Bending Anisotropy

The energetics of B-DNA bending toward the major and minor grooves were quantified by free energy simulations at four different KCl concentrations. Increased [KCl] led to more flexible DNA, with persistence lengths that agreed well with experimental values. At all salt concentrations, major groove bending was preferred, although preferences for major and minor groove bending were similar for the A-tract containing sequence. Since the phosphate repulsions and DNA internal energy favored minor groove bending, the preference for major groove bending was thought to originate from differences in solvation. Water in the minor groove was tighter bound than water in the major groove, and harder to displace than major groove water, which favored the compression of the major groove upon bending. Higher [KCl] decreased the persistence length for both major and minor groove bending but did not greatly affect the free energy spacing between the minor and major groove bending curves. For sequences without A-tracts, salt affected major and minor bending to nearly the same degree, and did not change the preference for major groove bending. For the A-tract containing sequence, an increase in salt concentration decreased the already small energetic difference between major and minor groove bending. Since salts did not significantly affect the relative differences in bending energetics and hydration, it is likely that the increased bending flexibilities upon salt increase are simply due to screening.

Consecutive non-natural PZ nucleobase pairs in DNA impact helical structure as seen in 50 μs molecular dynamics simulations

Z

Little is known about the influence of multiple consecutive 'non-standard' ( , 6-amino-5-nitro-2(1H)-pyridone, and , 2-amino-imidazo[1,2-a]-1,3,5-triazin-4(8H)-one) nucleobase pairs on the structural parameters of duplex DNA. nucleobase pairs follow standard rules for Watson-Crick base pairing but have rearranged hydrogen bonding donor and acceptor groups. Using the X-ray crystal structure as a starting point, we have modeled the motions of a DNA duplex built from a self-complementary oligonucleotide (5΄-CTTATPPPZZZATAAG-3΄) in water over a period of 50 μs and calculated DNA local parameters, step parameters, helix parameters, and major/minor groove widths to examine how the presence of multiple, consecutive nucleobase pairs might impact helical structure. In these simulations, the -containing DNA duplex exhibits a significantly wider major groove and greater average values of stagger, slide, rise, twist and h-rise than observed for a 'control' oligonucleotide in which nucleobase pairs are replaced by . The molecular origins of these structural changes are likely associated with at least two differences between and . First, the electrostatic properties of differ from in terms of density distribution and dipole moment. Second, differences are seen in the base stacking of pairs in dinucleotide steps, arising from energetically favorable stacking of the nitro group in with π-electrons of the adjacent base.

Revised RNA Dihedral Parameters for the Amber Force Field Improve RNA Molecular Dynamics

The backbone dihedral parameters of the Amber RNA force field were improved by fitting using multiple linear regression to potential energies determined by quantum chemistry calculations. Five backbone and four glycosidic dihedral parameters were fit simultaneously to reproduce the potential energies determined by a high-level density functional theory calculation (B97D3 functional with the AUG-CC-PVTZ basis set). Umbrella sampling was used to determine conformational free energies along the dihedral angles, and these better agree with the population of conformations observed in the protein data bank for the new parameters than for the conventional parameters. Molecular dynamics simulations performed on a set of hairpin loops, duplexes and tetramers with the new parameter set show improved modeling for the structures of tetramers CCCC, CAAU, and GACC, and an RNA internal loop of noncanonical pairs, as compared to the conventional parameters. For the tetramers, the new parameters largely avoid the incorrect intercalated structures that dominate the conformational samples from the conventional parameters. For the internal loop, the major conformation solved by NMR is stable with the new parameters, but not with the conventional parameters. The new force field performs similarly to the conventional parameters for the UUCG and GCAA hairpin loops and the [U(UA)6A]2 duplex.

Design and characterization of symmetric nucleic acids via molecular dynamics simulations

Asymmetry (5'→3') associated with each strand of the deoxyribonucleic acid (DNA) is inherent in the sugar-phosphate backbone connectivity and is essential for replication and transcription. We note that this asymmetry is due to one single chemical bond (C_3' to C_2' ) in each nucleotide unit, and the absence of this bond results in directionally symmetric nucleic acids. We also discovered that creation of an extra chemical bond (C_5' to C_2' ) can lead to a symmetric backbone. Keeping their potential synthetic and therapeutic interest in mind, we designed a few novel symmetric nucleic acids. We investigated their conformational stability and flexibility via detailed all atom explicit solvent 100-ns long molecular dynamics simulations and compared the resulting structures with that of regular B-DNA. Quite interestingly, some of the symmetric nucleic acids retain the overall double helical structure indicating their potential for integration in physiological DNA without causing major structural perturbations.

Influence of BII Backbone Substates on DNA Twist: A Unified View and Comparison of Simulation and Experiment for All 136 Distinct Tetranucleotide Sequences

Reliable representation of the B-DNA base-pair step twist is one of the crucial requirements for theoretical modeling of DNA supercoiling and other biologically relevant phenomena in B-DNA. It has long been suspected that the twist is inaccurately described by current empirical force fields. Unfortunately, comparison of simulation results with experiments is not straightforward because of the presence of BII backbone substates, whose populations may differ in experimental and simulation ensembles. In this work, we provide a comprehensive view of the effect of BII substates on the overall B-DNA helix twist and show how to reliably compare twist values from experiment and simulation in two scenarios. First, for longer DNA segments freely moving in solution, we show that sequence-averaged twists of different BI/BII ensembles can be compared directly because of approximate cancellation of the opposing BII effects. Second, for sequence-specific data, such as a particular base-pair step or tetranucleotide twist, can be compared only for a clearly defined BI/BII backbone conformation. For the purpose of force field testing, we designed a compact set of fourteen 22-base-pair B-DNA duplexes (Set 14) containing all 136 distinct tetranucleotide sequences and carried out a total of 84 μs of molecular dynamics simulations, primarily with the OL15 force field. Our results show that the ff99bsc0εζ_OL1χ_OL4, parmbsc1, and OL15 force fields model the B-DNA helical twist in good agreement with X-ray and minicircle ligation experiments. The comprehensive understanding obtained regarding the effect of BII substates on the base-pair step geometry should aid meaningful comparisons of various conformational ensembles in future research.

Anisotropy of B-DNA Groove Bending

DNA bending is critical for DNA packaging, recognition, and repair, and occurs toward either the major or the minor groove. The anisotropy of B-DNA groove bending was quantified for eight DNA sequences by free energy simulations employing a novel reaction coordinate. The simulations show that bending toward the major groove is preferred for non-A-tracts while the A-tract has a high tendency of bending toward the minor groove. Persistence lengths were generally larger for bending toward the minor groove, which is thought to originate from differences in groove hydration. While this difference in stiffness is one of the factors determining the overall preference of bending direction, the dominant contribution is shown to be a free energy offset between major and minor groove bending. The data suggests that, for the A-tract, this offset is largely determined by inherent structural properties, while differences in groove hydration play a large role for non-A-tracts. By quantifying the energetics of DNA groove bending and rationalizing the origins of the anisotropy, the calculations provide important new insights into a key biological process.

Assessing the Current State of Amber Force Field Modifications for DNA

The utility of molecular dynamics (MD) simulations to model biomolecular structure, dynamics, and interactions has witnessed enormous advances in recent years due to the availability of optimized MD software and access to significant computational power, including GPU multicore computing engines and other specialized hardware. This has led researchers to routinely extend conformational sampling times to the microsecond level and beyond. The extended sampling time has allowed the community not only to converge conformational ensembles through complete sampling but also to discover deficiencies and overcome problems with the force fields. Accuracy of the force fields is a key component, along with sampling, toward being able to generate accurate and stable structures of biopolymers. The Amber force field for nucleic acids has been used extensively since the 1990s, and multiple artifacts have been discovered, corrected, and reassessed by different research groups. We present a direct comparison of two of the most recent and state-of-the-art Amber force field modifications, bsc1 and OL15, that focus on accurate modeling of double-stranded DNA. After extensive MD simulations with five test cases and two different water models, we conclude that both modifications are a remarkable improvement over the previous bsc0 force field. Both force field modifications show better agreement when compared to experimental structures. To ensure convergence, the Drew–Dickerson dodecamer (DDD) system was simulated using 100 independent MD simulations, each extended to at least 10 μs, and the independent MD simulations were concatenated into a single 1 ms long trajectory for each combination of force field and water model. This is significantly beyond the time scale needed to converge the conformational ensemble of the internal portions of a DNA helix absent internal base pair opening. Considering all of the simulations discussed in the current work, the MD simulations performed to assess and validate the current force fields and water models aggregate over 14 ms of simulation time. The results suggest that both the bsc1 and OL15 force fields render average structures that deviate significantly less than 1 Å from the average experimental structures. This can be compared to similar but less exhaustive simulations with the CHARMM 36 force field that aggregate to the ∼90 μs time scale and also perform well but do not produce structures as close to the DDD NMR average structures (with root-mean-square deviations of 1.3 Å) as the newer Amber force fields. On the basis of these analyses, any future research involving double-stranded DNA simulations using the Amber force fields should employ the bsc1 or OL15 modification.

Transitions of Double-Stranded DNA Between the A- and B-Forms

The structure of double-stranded DNA (dsDNA) is sensitive to solvent conditions. In solution, B-DNA is the favored conformation under physiological conditions, while A-DNA is the form found under low water activity. The A-form is induced locally in some protein-DNA complexes, and repeated transitions between the B- and A-forms have been proposed to generate the forces used to drive dsDNA into viral capsids during genome packaging. Here, we report analyses on previous molecular dynamics (MD) simulations on B-DNA, along with new MD simulations on the transition from A-DNA to B-DNA in solution. We introduce the A-B Index (ABI), a new metric along the A-B continuum, to quantify our results. When A-DNA is placed in an equilibrated solution at physiological ionic strength, there is no energy barrier to the transition to the B-form, which begins within about 1 ns. The transition is essentially complete within 5 ns, although occasionally a stretch of a few base pairs will remain A-like for up to ∼10 ns. A comparison of four sequences with a range of predicted A-phobicities shows that more A-phobic sequences make the transition more rapidly than less A-phobic sequences. Simulations on dsDNA with a region of roughly one turn locked in the A-form allow us to characterize the A/B junction, which has an average bend angle of 20-30°. Fluctuations in this angle occur with characteristic times of about 10 ns.

Computational Approach to Explore the B/A Junction Free Energy in DNA

Protein-DNA interactions induce conformational changes in DNA such as B- to A-form transitions at a local level. Such transitions are associated with a junction free energy cost at the boundary of two different conformations in a DNA molecule. In this study, we performed umbrella sampling simulations to find the free energy values of the B-A transition at the dinucleotide and trinucleotide level of DNA. Using a combination of dinucleotide and trinucleotide free energy costs obtained from simulations, we calculated the B/A junction free energy. Our study shows that the B/A junction free energy is 0.52 kcal mol(-1) for the A-philic GG step and 1.59 kcal mol(-1) for the B-philic AA step. This observation is in agreement with experimentally derived values. After excluding junction effects, we obtained an absolute free energy cost for the B- to A-form conversion for all the dinucleotide steps. These absolute free energies may be used for predicting the propensity of structural transitions in DNA.