Mesoscopic model for free-energy-landscape analysis of DNA sequences

Bubble Relaxation Dynamics in Homopolymer DNA Sequences

Understanding the inherent timescales of large bubbles in DNA is critical to a thorough comprehension of its physicochemical characteristics, as well as their potential role on helix opening and biological function. In this work, we employ the coarse-grained Peyrard-Bishop-Dauxois model of DNA to study relaxation dynamics of large bubbles in homopolymer DNA, using simulations up to the microsecond time scale. By studying energy autocorrelation functions of relatively large bubbles inserted into thermalised DNA molecules, we extract characteristic relaxation times from the equilibration process for both adenine-thymine (AT) and guanine-cytosine (GC) homopolymers. Bubbles of different amplitudes and widths are investigated through extensive statistics and appropriate fittings of their relaxation. Characteristic relaxation times increase with bubble amplitude and width. We show that, within the model, relaxation times are two orders of magnitude longer in GC sequences than in AT sequences. Overall, our results confirm that large bubbles leave a lasting impact on the molecule's dynamics, for times between 0.5-500 ns depending on the homopolymer type and bubble shape, thus clearly affecting long-time evolutions of the molecule.

Nonlinear physics opens a new paradigm for accurate transcription start site prediction

There is evidence that DNA breathing (spontaneous opening of the DNA strands) plays a relevant role in the interactions of DNA with other molecules, and in particular in the transcription process. Therefore, having physical models that can predict these openings is of interest. However, this source of information has not been used before either in transcription start sites (TSSs) or promoter prediction. In this article, one such model is used as an additional information source that, when used by a machine learning (ML) model, improves the results of current methods for the prediction of TSSs. In addition, we provide evidence on the validity of the physical model, as it is able by itself to predict TSSs with high accuracy. This opens an exciting avenue of research at the intersection of statistical mechanics and ML, where ML models in bioinformatics can be improved using physical models of DNA as feature extractors.

Bubble lifetimes in DNA gene promoters and their mutations affecting transcription

Relative lifetimes of inherent double stranded DNA openings with lengths up to ten base pairs are presented for different gene promoters and corresponding mutants that either increase or decrease transcriptional activity in the framework of the Peyrard-Bishop-Dauxois model. Extensive microcanonical simulations are used with energies corresponding to physiological temperature. The bubble lifetime profiles along the DNA sequences demonstrate a significant reduction of the average lifetime at the mutation sites when the mutated promoter decreases transcription, while a corresponding enhancement of the bubble lifetime is observed in the case of mutations leading to increased transcription. The relative difference in bubble lifetimes between the mutated and wild type promoters at the position of mutation varies from 20% to more than 30% as the bubble length decreases.

Base pair fluctuations in helical models for nucleic acids

A statistical method is developed to estimate the maximum amplitude of the base pair fluctuations in a three dimensional mesoscopic model for nucleic acids. The base pair thermal vibrations around the helix diameter are viewed as a Brownian motion for a particle embedded in a stable helical structure. The probability to return to the initial position is computed, as a function of time, by integrating over the particle paths consistent with the physical properties of the model potential. The zero time condition for the first-passage probability defines the constraint to select the integral cutoff for various macroscopic helical conformations, obtained by tuning the twist, bending, and slide motion between adjacent base pairs along the molecule stack. Applying the method to a short homogeneous chain at room temperature, we obtain meaningful estimates for the maximum fluctuations in the twist conformation with ∼10.5 base pairs per helix turn, typical of double stranded DNA helices. Untwisting the double helix, the base pair fluctuations broaden and the integral cutoff increases. The cutoff is found to increase also in the presence of a sliding motion, which shortens the helix contour length, a situation peculiar of dsRNA molecules.

Effective Biophysical Modeling of Cell Free Transcription and Translation Processes

Transcription and translation are at the heart of metabolism and signal transduction. In this study, we developed an effective biophysical modeling approach to simulate transcription and translation processes. The model, composed of coupled ordinary differential equations, was tested by comparing simulations of two cell free synthetic circuits with experimental measurements generated in this study. First, we considered a simple circuit in which sigma factor 70 induced the expression of green fluorescent protein. This relatively simple case was then followed by a more complex negative feedback circuit in which two control genes were coupled to the expression of a third reporter gene, green fluorescent protein. Many of the model parameters were estimated from previous biophysical studies in the literature, while the remaining unknown model parameters for each circuit were estimated by minimizing the difference between model simulations and messenger RNA (mRNA) and protein measurements generated in this study. In particular, either parameter estimates from published studies were used directly, or characteristic values found in the literature were used to establish feasible ranges for the parameter estimation problem. In order to perform a detailed analysis of the influence of individual model parameters on the expression dynamics of each circuit, global sensitivity analysis was used. Taken together, the effective biophysical modeling approach captured the expression dynamics, including the transcription dynamics, for the two synthetic cell free circuits. While, we considered only two circuits here, this approach could potentially be extended to simulate other genetic circuits in both cell free and whole cell biomolecular applications as the equations governing the regulatory control functions are modular and easily modifiable. The model code, parameters, and analysis scripts are available for download under an MIT software license from the Varnerlab GitHub repository.

Distributions of bubble lifetimes and bubble lengths in DNA

We investigate the distribution of bubble lifetimes and bubble lengths in DNA at physiological temperature, by performing extensive molecular dynamics simulations with the Peyrard-Bishop-Dauxois (PBD) model, as well as an extended version (ePBD) having a sequence-dependent stacking interaction, emphasizing the effect of the sequences' guanine-cytosine (GC)/adenine-thymine (AT) content on these distributions. For both models we find that base pair-dependent (GC vs AT) thresholds for considering complementary nucleotides to be separated are able to reproduce the observed dependence of the melting temperature on the GC content of the DNA sequence. Using these thresholds for base pair openings, we obtain bubble lifetime distributions for bubbles of lengths up to ten base pairs as the GC content of the sequences is varied, which are accurately fitted with stretched exponential functions. We find that for both models the average bubble lifetime decreases with increasing either the bubble length or the GC content. In addition, the obtained bubble length distributions are also fitted by appropriate stretched exponential functions and our results show that short bubbles have similar likelihoods for any GC content, but longer ones are substantially more likely to occur in AT-rich sequences. We also show that the ePBD model permits more, longer-lived, bubbles than the PBD system.

Thermal versus mechanical unfolding in a model protein

Force spectroscopy techniques are often used to learn about the free energy landscape of single biomolecules, typically by recovering free energy quantities that, extrapolated to zero force, are compared to those measured in bulk experiments. However, it is not always clear how the information obtained from a mechanically perturbed system can be related to the information obtained using other denaturants since tensioned molecules unfold and refold along a reaction coordinate imposed by the force, which is not likely to be meaningful in its absence. Here, we explore this dichotomy by investigating the unfolding landscape of a model protein, which is unfolded first mechanically through typical force spectroscopy-like protocols and next thermally. When unfolded by nonequilibrium force extension and constant force protocols, we recover a simple two-barrier landscape as the protein reaches the extended conformation through a metastable intermediate. Interestingly, folding-unfolding equilibrium simulations at low forces suggested a totally different scenario, where this metastable state plays little role in the unfolding mechanism, and the protein unfolds through two competing pathways [R. Tapia-Rojo et al., J. Chem. Phys. 141, 135102 (2014)]. Finally, we use Markov state models to describe the configurational space of the unperturbed protein close to the critical temperature. The thermal dynamics is well understood by a one-dimensional landscape along an appropriate reaction coordinate, however it is very different from the mechanical picture. In this sense, the results of our protein model for the mechanical and thermal descriptions provide incompatible views of the folding/unfolding landscape of the system, and the estimated quantities to zero force result are hard to interpret.

Heterogeneity and chaos in the Peyrard-Bishop-Dauxois DNA model

We discuss the effect of heterogeneity on the chaotic properties of the Peyrard-Bishop-Dauxois nonlinear model of DNA. Results are presented for the maximum Lyapunov exponent and the deviation vector distribution. Different compositions of adenine-thymine (AT) and guanine-cytosine (GC) base pairs are examined for various energies up to the melting point of the corresponding sequence. We also consider the effect of the alternation index, which measures the heterogeneity of the DNA chain through the number of alternations between different types (AT or GC) of base pairs, on the chaotic behavior of the system. Biological gene promoter sequences have been also investigated, showing no distinct behavior of the maximum Lyapunov exponent.

Mesoscopic model and free energy landscape for protein-DNA binding sites: analysis of cyanobacterial promoters

The identification of protein binding sites in promoter sequences is a key problem to understand and control regulation in biochemistry and biotechnological processes. We use a computational method to analyze promoters from a given genome. Our approach is based on a physical model at the mesoscopic level of protein-DNA interaction based on the influence of DNA local conformation on the dynamics of a general particle along the chain. Following the proposed model, the joined dynamics of the protein particle and the DNA portion of interest, only characterized by its base pair sequence, is simulated. The simulation output is analyzed by generating and analyzing the Free Energy Landscape of the system. In order to prove the capacity of prediction of our computational method we have analyzed nine promoters of Anabaena PCC 7120. We are able to identify the transcription starting site of each of the promoters as the most populated macrostate in the dynamics. The developed procedure allows also to characterize promoter macrostates in terms of thermo-statistical magnitudes (free energy and entropy), with valuable biological implications. Our results agree with independent previous experimental results. Thus, our methods appear as a powerful complementary tool for identifying protein binding sites in promoter sequences.

Binding of specific proteins to particular sites in the DNA sequence is a fundamental issue for gene regulation in molecular biology and genetic engineering. A deep understanding of cell physiology requires the analysis of a plethora of genes involving characterization of their promoter architectures that determine their regulation and gene transcription. In order to locate the promoter elements of a given gene, experimental determination of its transcription start site (TSS) is required. This is an expensive, time-consuming task that, depending on our requirements, could be simplified using computational analysis as a first approach. Nevertheless, most computational methods lack a physical basis on the protein-DNA interaction mechanism. We adopt here this strategy, by using a simple model for protein-DNA interaction to find TSS in a bunch of cyanobacteria promoters. We make use of physical tools to characterize these TSS and to relate them with biological properties as the relative strength of the promoter. Our study shows how a model based on a coarse-grained description of a biomolecule can give valuable insight on its biological function.

Mesoscopic model parametrization of hydrogen bonds and stacking interactions of RNA from melting temperatures

Information about molecular interactions in DNA can be obtained from experimental melting temperature data by using mesoscopic statistical physics models. Here, we extend the technique to RNA and show that the new parameters correctly reproduce known properties such as the stronger hydrogen bonds of AU base pairs. We also were able to calculate a complete set of elastic constants for all 10 irreducible combinations of nearest neighbours (NNs). We believe that this is particularly useful as experimentally derived information about RNA elasticity is relatively scarce. The melting temperature prediction using the present model improves over those from traditional NN model, providing thus an alternative way to calculate these temperatures for RNA. Additionally, we calculated the site-dependent base pair oscillation to explain why RNA shows larger oscillation amplitudes despite having stronger AU hydrogen bonds.