On the absence of intrahelical DNA dynamics on the μs to ms timescale

Study of the Myosin Relay Helix Peptide by Molecular Dynamics Simulations, Pump-Probe and 2D Infrared Spectroscopy

The Davydov model was conjectured to describe how an amide I excitation created during ATP hydrolysis in myosin might be significant in providing energy to drive myosin's chemomechanical cycle. The free energy surfaces of the myosin relay helix peptide dissolved in 2,2,2-trifluoroethanol (TFE), determined by metadynamics simulations, demonstrate local minima differing in free energy by only ~2 kT, corresponding to broken and stabilized hydrogen bonds, respectively. Experimental pump-probe and 2D infrared spectroscopy were performed on the peptide dissolved in TFE. The relative heights of two peaks seen in the pump-probe data and the corresponding relative volumes of diagonal peaks seen in the 2D-IR spectra at time delays between 0.5 ps and 1 ps differ noticeably from what is seen at earlier or later time delays or in the linear spectrum, indicating that a vibrational excitation may influence the conformational state of this helix. Thus, it is possible that the presence of an amide I excitation may be a direct factor in the conformational state taken on by the myosin relay helix following ATP hydrolysis in myosin.

van der Waals Parameter Scanning with Amber Nucleic Acid Force Fields: Revisiting Means to Better Capture the RNA/DNA Structure through MD

Molecular dynamics simulations can be used in combination with experimental techniques to uncover the intricacies of biomolecular structure, dynamics, and the resulting interactions. However, many noncanonical nucleic acid structures have proven to be challenging to replicate in accurate agreement with experimental data, often attributed to known force field deficiencies. A common force field criticism is the handling of van der Waals (vdW) parameters, which have not been updated since the regular use of Ewald's methods became routine. This work dives into the effects of minute vdW radii shifts on RNA tetranucleotide, B-DNA, and Z-DNA model systems described by commonly used Amber force fields. Using multidimensional replica exchange molecular dynamics (M-REMD), the GACC RNA tetranucleotide demonstrated changes in the structural distribution between the NMR minor and anomalous structure populations based on the O2' vdW radii scanning. However, no significant change in the NMR Major conformation population was observed. There were minimal changes in the B-DNA structure but there were more substantial improvements in Z-DNA structural descriptions, specifically with the Tumuc1 force field. This occurred with both LJbb vdW radii adjustments and incorporation of the CUFIX nonbonded parameter modifications. Though the limited vdW modifications tested did not provide a universal fix to the challenge of simulating the various known nucleic acid structures, they do provide direction and a greater understanding for future force field development efforts.

Discovery of Novel Inhibitors of BRD4 for Treating Prostate Cancer: A Comprehensive Case Study for Considering Water Networks in Virtual Screening and Drug Design

Androgen receptor (AR) is the primary target for treating prostate cancer (PCa), which inevitably progresses due to drug-resistant mutations. Bromodomain-containing protein 4 (BRD4) has been a new potential drug target for PCa treatment. Herein, we report the rational design and discovery of novel BRD4 inhibitors through computer-aided drug design (CADD), and a hit compound SQ-1 (IC50 = 676 nM) was identified by structure-based virtual screening (SBVS) with the conserved water network. To optimize the structure of SQ-1, the free energy landscape was constructed, and the binding mechanism was explored by characterizing the water profile and the dissociation mechanism. Finally, the compound SQ-17 with improved inhibitory activity (IC50 < 100 nM) was discovered, which showed potent antiproliferative activity against LNCaP. These data highlighted a successful attempt to identify and optimize a small molecule by comprehensive CADD application and provided essential clues for developing novel therapeutics for PCa treatment.

Dynamic geometry design of cyclic peptide architectures for RNA structure

Designing inhibitors for RNA is still challenging due to the bottleneck of maintaining the binding interaction of inhibitor-RNA accompanied by subtle RNA flexibility. Thus, the current approach usually needs to screen thousands of candidate inhibitors for binding. Here, we propose a dynamic geometry design approach to enrich the hits with only a tiny pool of designed geometrically compatible scaffold candidates. First, our method uses graph-based tree decomposition to explore the complementarity rigid binding cyclic peptide and design the amino acid side chain length and charge to fit the RNA pocket. Then, we perform an energy-based dynamical network algorithm to optimize the inhibitor-RNA hydrogen bonds. Dynamic geometry-guided design yields successful inhibitors with low micromolar binding affinity scaffolds and experimentally competes with the natural RNA chaperone. The results indicate that the dynamic geometry method yields higher efficiency and accuracy than traditional methods. The strategy could be further optimized to design the length and chirality by adopting nonstandard amino acids and facilitating RNA engineering for biological or medical applications.

Mechanical codes of chemical-scale specificity in DNA motifs

In gene transcription, certain sequences of double-stranded (ds)DNA play a vital role in nucleosome positioning and expression initiation. That dsDNA is deformed to various extents in these processes leads us to ask: Could the genomic DNA also have sequence specificity in its chemical-scale mechanical properties? We approach this question using statistical machine learning to determine the rigidity between DNA chemical moieties. What emerges for the polyA, polyG, TpA, and CpG sequences studied here is a unique trigram that contains the quantitative mechanical strengths between bases and along the backbone. In a way, such a sequence-dependent trigram could be viewed as a DNA mechanical code. Interestingly, we discover a compensatory competition between the axial base-stacking interaction and the transverse base-pairing interaction, and such a reciprocal relationship constitutes the most discriminating feature of the mechanical code. Our results also provide chemical-scale understanding for experimental observables. For example, the long polyA persistence length is shown to have strong base stacking while its complement (polyAc) exhibits high backbone rigidity. The mechanical code concept enables a direct reading of the physical interactions encoded in the sequence which, with further development, is expected to shed new light on DNA allostery and DNA-binding drugs.

In silico prediction of mutation sites for anthranilate synthase from Serratia marcesens to deregulate tryptophan feedback inhibition

Allosteric feedback inhibition of the committed step in amino acid biosynthetic pathways is a major concern for production of amino acids at industrial scale. Anthranilate synthase (AS) catalyzes the first reaction of tryptophan biosynthetic pathway found in microorganisms and is feedback inhibited by its own product i.e. tryptophan. Here, we identified new mutant sites in AS using computational mutagenesis approach. MD simulations (20 ns) followed by MMPBSA and per residue decomposition energy analysis identified seven amino acid residues with best binding affinity for tryptophan. All 19 mutant structures were generated for each identified amino acid residue followed by simulation to evaluate effect of mutation on protein stability. Later, molecular docking studies were employed to generate mutant-tryptophan complex and structures with binding energies (kcal/mol) much higher than wild-type AS were selected. Finally, two mutants i.e., S37W and S37H were identified on the basis of positive binding scores and loss of tryptophan binding inside pocket. Further, MD simulations run for 200 ns were performed over these mutant-tryptophan complexes followed by RMSD, RMSF, radius of gyration , solvent accessible surface area , intra-protein hydrogen bond numbers, principal component analysis, free energy landscape (FEL) and secondary structure analysis to rationale effect of mutations on stability of protein. Cross correlation analysis of mutant site amino acids (S37W) with key residues of catalytic site (G325, T326, H395 and G482) was done to evaluate the effect of mutations on catalytic site conformation. Current computational mutagenesis approach predicted two mutants S37W and S37H with proposed deregulated feedback inhibition by tryptophan and retained catalytic activity.Communicated by Ramaswamy H. Sarma.

Computational Modeling of DNA 3D Structures: From Dynamics and Mechanics to Folding

DNA carries the genetic information required for the synthesis of RNA and proteins and plays an important role in many processes of biological development. Understanding the three-dimensional (3D) structures and dynamics of DNA is crucial for understanding their biological functions and guiding the development of novel materials. In this review, we discuss the recent advancements in computer methods for studying DNA 3D structures. This includes molecular dynamics simulations to analyze DNA dynamics, flexibility, and ion binding. We also explore various coarse-grained models used for DNA structure prediction or folding, along with fragment assembly methods for constructing DNA 3D structures. Furthermore, we also discuss the advantages and disadvantages of these methods and highlight their differences.

Selectivity and Ranking of Tight-Binding JAK-STAT Inhibitors Using Markovian Milestoning with Voronoi Tessellations

Janus kinases (JAK), a group of proteins in the nonreceptor tyrosine kinase (NRTKs) family, play a crucial role in growth, survival, and angiogenesis. They are activated by cytokines through the Janus kinase-signal transducer and activator of a transcription (JAK-STAT) signaling pathway. JAK-STAT signaling pathways have significant roles in the regulation of cell division, apoptosis, and immunity. Identification of the V617F mutation in the Janus homology 2 (JH2) domain of JAK2 leading to myeloproliferative disorders has stimulated great interest in the drug discovery community to develop JAK2-specific inhibitors. However, such inhibitors should be selective toward JAK2 over other JAKs and display an extended residence time. Recently, novel JAK2/STAT5 axis inhibitors (N-(1H-pyrazol-3-yl)pyrimidin-2-amino derivatives) have displayed extended residence times (hours or longer) on target and adequate selectivity excluding JAK3. To facilitate a deeper understanding of the kinase-inhibitor interactions and advance the development of such inhibitors, we utilize a multiscale Markovian milestoning with Voronoi tessellations (MMVT) approach within the Simulation-Enabled Estimation of Kinetic Rates v.2 (SEEKR2) program to rank order these inhibitors based on their kinetic properties and further explain the selectivity of JAK2 inhibitors over JAK3. Our approach investigates the kinetic and thermodynamic properties of JAK-inhibitor complexes in a user-friendly, fast, efficient, and accurate manner compared to other brute force and hybrid-enhanced sampling approaches.

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.

Flexibility of flanking DNA is a key determinant of transcription factor affinity for the core motif

Selective gene regulation is mediated by recognition of specific DNA sequences by transcription factors (TFs). The extremely challenging task of searching out specific cognate DNA binding sites among several million putative sites within the eukaryotic genome is achieved by complex molecular recognition mechanisms. Elements of this recognition code include the core binding sequence, the flanking sequence context, and the shape and conformational flexibility of the composite binding site. To unravel the extent to which DNA flexibility modulates TF binding, in this study, we employed experimentally guided molecular dynamics simulations of ternary complex of closely related Hox heterodimers Exd-Ubx and Exd-Scr with DNA. Results demonstrate that flexibility signatures embedded in the flanking sequences impact TF binding at the cognate binding site. A DNA sequence has intrinsic shape and flexibility features. While shape features are localized, our analyses reveal that flexibility features of the flanking sequences percolate several basepairs and allosterically modulate TF binding at the core. We also show that lack of flexibility in the motif context can render the cognate site resistant to protein-induced shape changes and subsequently lower TF binding affinity. Overall, this study suggests that flexibility-guided DNA shape, and not merely the static shape, is a key unexplored component of the complex DNA-TF recognition code.

Insights into the binding mode of AS1411 aptamer to nucleolin

AS1411 aptamer can function as a recognition probe to detect the cell surface nucleolin overexpressed in cancer cells, however, little is known about their binding process. This study proposed a feasible binding mode for the first time and provided atomic-level descriptions for the high affinity and specific binding of AS1411. The binding pose predicted by docking was screened using knowledge-based criteria, and a microsecond molecular dynamics (MD) simulation showed the stable existence of the predicted structure in the solution. Structural analysis shows that the unique capping of the 5′ end of AS1411 provides the specific binding with RBD1, and the interactions of hydrogen bond, salt bridge, and water-mediated network between AS1411 and RBD1,2 stabilize the binding. The calculation of per-residue decomposition emphasizes the dominant contribution of van der Waals energy and critical residues are screened. Our study provides the molecular basis of this specific binding and can guide rational AS1411-based aptamers design. Further insights require tight collaborations between the experiments and in silico studies.

Transient Hoogsteen Base Pairs Observed in Unbiased Molecular Dynamics Simulations of DNA

Duplex DNA is modeled as canonical B-DNA displaying the characteristic Watson-Crick base pairs. A less common and short-lived pairing of the nucleobases is the Hoogsteen (HG) conformation. The low population of the HG base pairs (<1%) necessitates extended sampling times in order to analyze through unbiased molecular dynamics (MD) simulations. We have discovered that with extended sampling times using multiple independent copies of an 18-mer sequence, the MD trajectories reproduce the expected and transient HG base pairing. Consistent with experiment, the percentage of the HG events are within the range of ∼0.1-1.0% over the combined aggregate sampling time of more than 3.6 ms. We present the reliability of the current AMBER set of nucleic acid force fields and tools to accurately simulate naturally occurring base-pairing and opening events without any bias or restraints. The mechanism consists of base pair fraying, flipping of the purine, and reformation with HG base pairs.

1H NMR Chemical Exchange Techniques Reveal Local and Global Effects of Oxidized Cytosine Derivatives

5-Carboxycytosine (5caC) is a rare epigenetic modification found in nucleic acids of all domains of life. Despite its sparse genomic abundance, 5caC is presumed to play essential regulatory roles in transcription, maintenance and base-excision processes in DNA. In this work, we utilize nuclear magnetic resonance (NMR) spectroscopy to address the effects of 5caC incorporation into canonical DNA strands at multiple pH and temperature conditions. Our results demonstrate that 5caC has a pH-dependent global destabilizing and a base-pair mobility enhancing local impact on dsDNA, albeit without any detectable influence on the ground-state B-DNA structure. Measurement of hybridization thermodynamics and kinetics of 5caC-bearing DNA duplexes highlighted how acidic environment (pH 5.8 and 4.7) destabilizes the double-stranded structure by ∼10–20 kJ mol^–1 at 37 °C when compared to the same sample at neutral pH. Protonation of 5caC results in a lower activation energy for the dissociation process and a higher barrier for annealing. Studies on conformational exchange on the microsecond time scale regime revealed a sharply localized base-pair motion involving exclusively the modified site and its immediate surroundings. By direct comparison with canonical and 5-formylcytosine (5fC)-edited strands, we were able to address the impact of the two most oxidized naturally occurring cytosine derivatives in the genome. These insights on 5caC’s subtle sensitivity to acidic pH contribute to the long-standing questions of its capacity as a substrate in base excision repair processes and its purpose as an independent, stable epigenetic mark.

Searching for Low Probability Opening Events in a DNA Sliding Clamp

The β subunit of E. coli DNA polymererase III is a DNA sliding clamp associated with increasing the processivity of DNA synthesis. In its free form, it is a circular homodimer structure that can accomodate double-stranded DNA in a nonspecific manner. An open state of the clamp must be accessible before loading the DNA. The opening mechanism is still a matter of debate, as is the effect of bound DNA on opening/closing kinetics. We use a combination of atomistic, coarse-grained, and enhanced sampling strategies in both explicit and implicit solvents to identify opening events in the sliding clamp. Such simulations of large nucleic acid and their complexes are becoming available and are being driven by improvements in force fields and the creation of faster computers. Different models support alternative opening mechanisms, either through an in-plane or out-of-plane opening event. We further note some of the current limitations, despite advances, in modeling these highly charged systems with implicit solvent.

Structural basis for molecular recognition of G protein-coupled estrogen receptor by selected bisphenols

Complicated ligand-dependent signaling pathways of bisphenol A (BPA) and its analogues involve not only intranuclear estrogen receptor but also membrane receptor G protein-coupled estrogen receptor (GPER). However, the structural basis for molecular recognition of GPER by the environmental chemicals remains unknown. To reveal the structural dependence of GPER recognition by bisphenols, a systematic molecular dynamics simulation study was performed for selected bisphenols with different electron hybrid orbitals and substituents on their C atoms connecting two phenol rings. BPA was used as a control, bisphenol C(BPC) as an example for a connecting C with sp² hybrid orbitals to provide more ligand rigidity, bisphenol E(BPE) and bisphenol F(BPF) for decreased steric hindrance and hydrophobicity around the connecting C, and bisphenol B(BPB) and bisphenol AF(BPAF) for increased hydrophobicity and steric hindrance. All the tested bisphenols can bind with GPER at its classic orthosteric site to obtain GPER-ligand complexes, while van der Waals interactions and direct inter-molecular electrostatic energies provide the driving forces for ligand binding. Bulky substituents and structural rigidity of the connecting C dramatically impair hydrogen bonding between GPER and the bisphenols, which results in decreased contribution of both favorable intermolecular hydrogen bonds and unfavorable polar solvation effect to complex stability of BPB and BPC since decreased number of key residues is expected. Increase in substituent lipophilicity enhances the van der Waals interactions and favorable non-polar solvation effect. The six bisphenols of high structural similarity shared two key recognition residues, Leu137^TM3 and Trp272^TM6, the latter of which was in the highly conserved CWxP motif of TM6 and has been reported as key residue for G protein-coupled receptor activation. Based on the obtained knowledge, GPER affinity and relevant toxicity of BPA alternatives can be easily predicted, and the calculated binding free energies are consistent with the available experimental observations.

Recombinant laccase rPOXA 1B real-time, accelerated and molecular dynamics stability study

Background

Laccases (EC 1.10.3.2) are multi-copper oxidoreductases with great biotechnological importance due to their high oxidative potential and utility for removing synthetic dyes, oxidizing phenolic compounds, and degrading pesticides, among others.

Methods

A real-time stability study (RTS) was conducted for a year, by using enzyme concentrates from 3 batches (L1, L3, and L4). For which, five temperatures 243.15, 277.15, 298.15, 303.15, 308.15, and 313.15 K were assayed. Using RTS data and the Arrhenius equation, we calculated the rPOXA 1B accelerated stability (AS). Molecular dynamics (MD) computational study results were very close to those obtained experimentally at four different temperatures 241, 278, 298, and 314 K.

Results

In the RTS, 101.16, 115.81, 75.23, 46.09, 5.81, and 4.83% of the relative enzyme activity were recovered, at respective assayed temperatures. AS study, showed that rPOXA 1B is stable at 240.98 ± 5.38, 277.40 ± 1.32 or 297.53 ± 3.88 K; with t_1/2 values of 230.8, 46.2, and 12.6 months, respectively. Kinetic and thermodynamic parameters supported the high stability of rPOXA 1B, with an E_d value of 41.40 KJ mol^− 1, a low variation of K_M and V_max, at 240.98 ± 5.38, and 297.53 ± 3.88 K, and ∆G values showing deactivation reaction does not occur. The MD indicates that fluctuations in loop, coils or loops with hydrophilic or intermediate polarity amino acids as well as in some residues of POXA 1B 3D structure, increases with temperature; changing from three fluctuating residues at 278 K to six residues at 298 K, and nine residues at 314 K.

Conclusions

Laccase rPOXA 1B demonstrated experimentally and computationally to be a stable enzyme, with t_1/2 of 230.8, 46.2 or 12.6 months, if it is preserved impure without preservatives at temperatures of 240.98 ± 5.38, 277.40 ± 1.32 or 297.53 ± 3.88 K respectively; this study could be of great utility for large scale producers.

Supplementary Information

The online version contains supplementary material available at 10.1186/s12896-021-00698-3.

Ethidium bromide interactions with DNA: an exploration of a classic DNA-ligand complex with unbiased molecular dynamics simulations

Visualization of double stranded DNA in gels with the binding of the fluorescent dye ethidium bromide has been a basic experimental technique in any molecular biology laboratory for >40 years. The interaction between ethidium and double stranded DNA has been observed to be an intercalation between base pairs with strong experimental evidence. This presents a unique opportunity for computational chemistry and biomolecular simulation techniques to benchmark and assess their models in order to see if the theory can reproduce experiments and ultimately provide new insights. We present molecular dynamics simulations of the interaction of ethidium with two different double stranded DNA models. The first model system is the classic sequence d(CGCGAATTCGCG)₂ also known as the Drew–Dickerson dodecamer. We found that the ethidium ligand binds mainly stacked on, or intercalated between, the terminal base pairs of the DNA with little to no interaction with the inner base pairs. As the intercalation at the terminal CpG steps is relatively rapid, the resultant DNA unwinding, rigidification, and increased stability of the internal base pair steps inhibits further intercalation. In order to reduce these interactions and to provide a larger groove space, a second 18-mer DNA duplex system with the sequence d(GCATGAACGAACGAACGC) was tested. We computed molecular dynamics simulations for 20 independent replicas with this sequence, each with ∼27 μs of sampling time. Results show several spontaneous intercalation and base-pair eversion events that are consistent with experimental observations. The present work suggests that extended MD simulations with modern DNA force fields and optimized simulation codes are allowing the ability to reproduce unbiased intercalation events that we were not able to previously reach due to limits in computing power and the lack of extensively tested force fields and analysis tools.

Dual-Site Binding of Quaternary Ammonium Ions as Internal K+-Ion Channel Blockers: Nonclassical (C-H···O) H Bonding vs Dispersive (C-H···H-C) Interaction

A molecular-level study of the influence of the alkyl chain length of quaternary ammonium ions (QAs) on the blocking action and the mode of binding with the bacterial KcsA K⁺-ion channel is carried out by molecular dynamics (MD) simulations as well as quantum mechanics/molecular mechanics (QM/MM) methods. The present work unveils distinct modes of binding for different QAs, due to differences in size and hydrophobicity. The QAs bind near the channel gate as well as at the central cavity, leading to a possible dual-site blocking action. Small-sized tetraethylammonium (TEA) and tetrabutylammonium (TBA) ions enter inside the channel cavity in the open state of KcsA but bind strongly in the closed state. TEA binds to the polar hydroxyl group of threonine residues situated at the channel gate via nonclassical H-bonding interaction (C-H···O), while TBA binds to a second binding site, the central cavity, with hydrophobic benzyl and sec-butyl side chains of phenylalanine and isoleucine residues via alkyl-π and hydrophobic interactions (C-H···H-C). On the contrary, large tetrahexylammonium (THA) and tetraoctylammonium (TOA) ions bind the hydrophobic side-chain methyl and isopropyl of alanine and valine at the channel gate both in the open and closed states, thereby restricting the free movement of large QAs toward the center of the cavity. However, the binding to the hydrophobic benzyl and sec-butyl side chains of phenylalanine and isoleucine residues in the closed state is thermodynamically preferable. Also, the binding energy is found to increase with an increase in the alkyl chain length from ethyl (-16.4 kcal/mol) to octyl (-65.5 kcal/mol), due to an almost linear increase in dispersive interaction.

Ancillary Ligand in Ternary CuII Complexes Guides Binding Selectivity toward Minor-Groove DNA

Copper-containing compounds known as Casiopeínas are biologically active molecules which show promising antineoplastic effects against several cancer types. Two possible hypotheses regarding the mode of action of the Casiopeínas have emerged from the experimental evidence: the generation of reactive oxygen species or the ability of the compounds to bind and interact with nucleic acids. Using robust molecular dynamics simulations, we investigate the interaction of four different Casiopeínas with the DNA duplex d(GCACGAACGAACGAACGC). The studied copper complexes contain either 4–7- or 5–6-substituted dimethyl phenanthroline as the primary ligand and either glycinate or acetylacetonate as the secondary ligand. For statistical significance and to reduce bias in the simulations, four molecules of each copper compound were manually placed at a distance of 10 Å away from the DNA and 20 independent molecular dynamics simulations were performed, each reaching at least 30 μs. This time scale allows us to reproduce expected DNA terminal base-pair fraying and also to observe intercalation/base-pair eversion events generated by the compounds interacting with DNA. The results reveal that the secondary ligand is the guide toward the mode of binding between the copper complex and DNA in which glycinate prefers minor-groove binding and acetylacetonate produces base-pair eversion and intercalation. The Cu^II complexes containing glycinate interact within the DNA minor groove which are stabilized principally by the hydrogen bonds formed between the amino group of the aminoacidate moiety, whereas the compounds with the acetylacetonate do not present a stable network of hydrogen bonds and the ligand interactions enhance DNA breathing dynamics that result in base-pair eversion.

Screening of Natural Products Targeting SARS-CoV-2-ACE2 Receptor Interface - A MixMD Based HTVS Pipeline

The COVID-19 pandemic, caused by novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), is a severe global health crisis now. SARS-CoV-2 utilizes its Spike protein receptor-binding domain (S-protein) to invade human cell through binding to Angiotensin-Converting Enzyme 2 receptor (ACE2). S-protein is the key target for many therapeutics and vaccines. Potential S-protein-ACE2 fusion inhibitor is expected to block the virus entry into the host cell. In many countries, traditional practices, based on natural products (NPs) have been in use to slow down COVID-19 infection. In this study, a protocol was applied that combines mixed solvent molecular dynamics simulations (MixMD) with high-throughput virtual screening (HTVS) to search NPs to block SARS-CoV-2 entry into the human cell. MixMD simulations were employed to discover the most promising stable binding conformations of drug-like probes in the S-protein-ACE2 interface. Detected stable sites were used for HTVs of 612093 NPs to identify molecules that could interfere with the S-protein-ACE2 interaction. In total, 19 NPs were selected with rescoring model. These top-ranked NP-S-protein complexes were subjected to classical MD simulations for 300 ns (3 replicates of 100 ns) to estimate the stability and affinity of binding. Three compounds, ZINC000002128789, ZINC000002159944 and SN00059335, showed better stability in all MD runs, of which ZINC000002128789 was predicted to have the highest binding affinity, suggesting that it could be effective modulator in RBD-ACE2 interface to prevent SARS-CoV-2 infection. Our results support that NPs may provide tools to fight COVID-19.

Molecular dynamics simulation and QM/MM calculation reveal the selectivity mechanism of type I 1/2 kinase inhibitors: the effect of intramolecular H-bonds and conformational restriction for improved selectivity

Understanding the selectivity mechanisms of inhibitors towards highly similar proteins is extremely important work on the way to a new drug. Here, we aim to reveal the selectivity mechanisms of type I 1/2 kinase inhibitors towards p21-activated kinase (PAK4) and mitogen-activated protein kinase kinase kinase 14 (MAP3K14, NIK). PAK4, belonging to the serine/threonine protein kinases, is involved in cell signaling pathways and controls cellular functions and has received attention as an attractive drug target. The high sequence identity between PAK4 and NIK makes it challenging to design selective PAK4 inhibitors. In this work, computational methods including protein comparison, molecular docking, QM/MM, molecular dynamics simulations, and density functional theory (DFT) calculation were employed to explore the binding mechanisms of selective inhibitors against NIK and PAK4. The simulation results revealed the crucial factors accounting for selective inhibition of PAK4 over NIK, including different protein-ligand interactions, the positions and conformations of key residues, and the ligands flexibilities. This study will shed light on understanding the selectivity mechanisms of PAK4 and NIK inhibitors.

Asymmetric base-pair opening drives helicase unwinding dynamics

The opening of a Watson-Crick double helix is required for crucial cellular processes, including replication, repair, and transcription. It has long been assumed that RNA or DNA base pairs are broken by the concerted symmetric movement of complementary nucleobases. By analyzing thousands of base-pair opening and closing events from molecular simulations, here, we uncover a systematic stepwise process driven by the asymmetric flipping-out probability of paired nucleobases. We demonstrate experimentally that such asymmetry strongly biases the unwinding efficiency of DNA helicases toward substrates that bear highly dynamic nucleobases, such as pyrimidines, on the displaced strand. Duplex substrates with identical thermodynamic stability are thus shown to be more easily unwound from one side than the other, in a quantifiable and predictable manner. Our results indicate a possible layer of gene regulation coded in the direction-dependent unwindability of the double helix.

Lessons learned in atomistic simulation of double-stranded DNA: Solvation and salt concerns [Article v1.0]

Nucleic acids are highly charged macromolecules sensitive to their surroundings of water, salt, and other biomolecules. Molecular dynamics simulations with accurate biomolecular force fields provide a detailed atomistic view into DNA and RNA that has been useful to study the structure and dynamics of these molecules and their biological relevance. In this work we study the Drew-Dickerson dodecamer duplex with the sequence d(GCGCAATTGCGC)₂ in three different salt concentrations and using different monvalent salt types to detect possible structural influence. Overall, the DNA shows no major structural changes regardless of amount or type of monovalent ions used. Our results show that only at very high salt conditions (5M) is a small structural effect observed in the DNA duplex, which mainly consist of narrowing of the grooves due to increased residence of ions. We also present the importance of sampling time to achieve a converged ensemble, which is of major relevance in any simulation to avoid biased or non-meaningful results.

Sulfur-substitution-induced base flipping in the DNA duplex

Base flipping is widely observed in a number of important biological processes. The genetic codes deposited inside the DNA duplex become accessible to external agents upon base flipping. The sulfur substitution of guanine leads to thioguanine, which alters the thermodynamic stability of the GC base pairs and the GT mismatches. Experimental studies conclude that the sulfur substitution decreases the lifetime of the GC base pair. In this work, under three AMBER force fields for nucleotide systems, we firstly performed equilibrium and nonequilibrium free energy simulations to investigate the variation of the thermodynamic profiles in base flipping upon sulfur substitution. It is found that the bsc0 modification, the bsc1 modification and the OL15 modification of AMBER force fields are able to qualitatively describe the sulfur-substitution dependent behavior of the thermodynamics. However, only the two last-generation AMBER force fields are able to provide quantitatively correct predictions. The second computational study on the sulfur substitutions focused on the relative stability of the S6G-C base pair and the S6G-T mismatch. Two conflicting experimental observations were reported by the same authors. One suggested that the S6G-C base pair was more stable, while the other concludes that the S6G-T mismatch was more stable. We answered this question by constructing the free energy profiles along the base flipping pathway computationally.

Exploring potentially alternative non-canonical DNA duplex structures through simulation

Hopkins proposed an alternative and chirally distinct family of double-stranded DNA (dsDNA) models that have antiparallel chains with 5'→3' senses opposite to those of the right-handed Watson-Crick (WC) family. Termed configuration II, this family of dsDNA models contains both right-handed (II-R) and left-handed (II-L) forms, with Z-DNA as an example of the latter. Relative interstrand binding energies for six DNA duplex models, two each of configuration I-R (standard WC canonical B-DNA), II-R, and II-L for the duplex d(CGCGAATTCGCG), have been estimated under identical conditions using MM-PBSA analysis from molecular dynamics trajectories using three different AMBER force fields. These simulations support the stereo chemical soundness of configuration II dsDNA forms. Recent force fields (Barcelona Supercomputing Center [BSC] [bsc1] and Olomouc 2015 [OL15]) successfully render stable II-L structures, whereas the previous force field, bsc0, generated stable II-R structures, although with an energy difference between II-R and II-L of ∼30 kcal/mol. Communicated by Ramaswamy H. Sarma.

Energetics of base flipping at a DNA mismatch site confined at the latch constriction of α-hemolysin

Unique, two-state modulating current signatures are observed when a cytosine-cytosine mismatch pair is confined at the 2.4 nm latch constriction of the α-hemolysin (αHL) nanopore. We have previously speculated that the modulation is due to base flipping at the mismatch site. Base flipping is a biologically significant mechanism in which a single base is rotated out of the DNA helical stack by 180°. It is the mechanism by which enzymes are able to access bases for repair operations without disturbing the global structure of the helix. Here, temperature dependent ion channel recordings of individual double-stranded DNA duplexes inside αHL are used to derive thermodynamic (ΔH, ΔS) and kinetic (E_A) parameters for base flipping of a cytosine at an unstable cytosine-cytosine mismatch site. The measured activation energy for flipping a cytosine located at the latch of αHL out of the helix (18 ± 1 kcal mol^-1) is comparable to that previously reported for base flipping at mismatch sites from NMR measurements and potential mean force calculations. We propose that the αHL nanopore is a useful tool for measuring conformational changes in dsDNA at the single molecule level.

Molecular Modeling Applied to Nucleic Acid-Based Molecule Development

Molecular modeling by means of docking and molecular dynamics (MD) has become an integral part of early drug discovery projects, enabling the screening and enrichment of large libraries of small molecules. In the past decades, special emphasis was drawn to nucleic acid (NA)-based molecules in the fields of therapy, diagnosis, and drug delivery. Research has increased dramatically with the advent of the SELEX (systematic evolution of ligands by exponential enrichment) technique, which results in single-stranded DNA or RNA sequences that bind with high affinity and specificity to their targets. Herein, we discuss the role and contribution of docking and MD to the development and optimization of new nucleic acid-based molecules. This review focuses on the different approaches currently available for molecular modeling applied to NA interaction with proteins. We discuss topics ranging from structure prediction to docking and MD, highlighting their main advantages and limitations and the influence of flexibility on their calculations.

Ionic liquids make DNA rigid

Persistence length of double-stranded DNA (dsDNA) is known to decrease with an increase in ionic concentration of the solution. In contrast to this, here we show that the persistence length of dsDNA increases dramatically as a function of ionic liquid (IL) concentration. Using all atom explicit solvent molecular dynamics simulations and theoretical models, we present, for the first time, a systematic study to determine the mechanical properties of dsDNA in various hydrated ILs at different concentrations. We find that dsDNA in 50 wt % ILs have lower persistence length and stretch modulus in comparison to 80 wt % ILs. We further observe that both the persistence length and stretch modulus of dsDNA increase as we increase the concentration of ILs. The present trend of the stretch modulus and persistence length of dsDNA with IL concentration supports the predictions of the macroscopic elastic theory, in contrast to the behavior exhibited by dsDNA in monovalent salt. Our study further suggests the preferable ILs that can be used for maintaining DNA stability during long-term storage.

Dynamics of physiologically relevant noncanonical DNA structures: an overview from experimental and theoretical studies

DNA is a complex molecule with phenomenal inherent plasticity and the ability to form different hydrogen bonding patterns of varying stabilities. These properties enable DNA to attain a variety of structural and conformational polymorphic forms. Structurally, DNA can exist in single-stranded form or as higher-order structures, which include the canonical double helix as well as the noncanonical duplex, triplex and quadruplex species. Each of these structural forms in turn encompasses an ensemble of dynamically heterogeneous conformers depending on the sequence composition and environmental context. In vivo, the widely populated canonical B-DNA attains these noncanonical polymorphs during important cellular processes. While several investigations have focused on the structure of these noncanonical DNA, studying their dynamics has remained nontrivial. Here, we outline findings from some recent advanced experimental and molecular simulation techniques that have significantly contributed toward understanding the complex dynamics of physiologically relevant noncanonical forms of DNA.

AMOEBA Polarizable Atomic Multipole Force Field for Nucleic Acids

The AMOEBA polarizable atomic multipole force field for nucleic acids is presented. Valence and electrostatic parameters were determined from high-level quantum mechanical data, including structures, conformational energy, and electrostatic potentials, of nucleotide model compounds. Previously derived parameters for the phosphate group and nucleobases were incorporated. A total of over 35 μs of condensed-phase molecular dynamics simulations of DNA and RNA molecules in aqueous solution and crystal lattice were performed to validate and refine the force field. The solution and/or crystal structures of DNA B-form duplexes, RNA duplexes, and hairpins were captured with an average root-mean-squared deviation from NMR structures below or around 2.0 Å. Structural details, such as base pairing and stacking, sugar puckering, backbone and χ-torsion angles, groove geometries, and crystal packing interfaces, agreed well with NMR and/or X-ray. The interconversion between A- and B-form DNAs was observed in ethanol-water mixtures at 328 K. Crystal lattices of B- and Z-form DNA and A-form RNA were examined with simulations. For the RNA tetraloop, single strand tetramers, and HIV TAR with 29 residues, the simulated conformational states, 3 J-coupling, nuclear Overhauser effect, and residual dipolar coupling data were compared with NMR results. Starting from a totally unstacked/unfolding state, the rCAAU tetranucleotide was folded into A-form-like structures during ∼1 μs molecular dynamics simulations.

Structural and dynamical instability of DNA caused by high occurrence of d5SICS and dNaM unnatural nucleotides

In vivo inclusion of unnatural base pairs (UBPs) into functional DNA was recently reported for compounds 2,6-dimethyl-2H-isoquiniline-1-thione (d5SICS, X) and 2-methoxy-3-methylnaphthalene (dNaM, Y) in a modified E. coli strand, for which high-fidelity replication was observed. Yet, little is known about possible genetic consequences they have in largely substituted DNA. Using a converged microsecond long molecular dynamics (MD) simulation of the sequences 5'-GCGCAAXTTGCGC-3' and 5'-GCGCXAXTXGCGC-3', where X represents the UBP, we show that in the system with only a single XY UBP pair present, the global RMS deviation from canonical B-DNA in our control simulations is ∼3 Å and a fully converged ensemble is achieved within 2 µs. With three UBPs, deviation increases to ∼5 Å and convergence is not achieved within 10 µs of sampling time. With five UBPs, no convergence is observed and the double helix collapses into a globular structure. A fully optimized structure of the trimer d(GXC) was obtained using the density functional theory method B97D/cc-pVTZ and resulted in an RMSD value of ∼2 Å when compared to the most populated structure obtained from the MD simulations. Their calculated interaction energy is -3.7 kcal mol^-1. It is thus unlikely that d5SICS and dNaM could be useful as tools in DNA manipulation. This theoretical methodology can be used in the de novo design of UBPs.

All-atom molecular dynamics simulations of spin labelled double and single-strand DNA for EPR studies

Prediction of motional EPR spectra of spin labelled DNA structures from fully atomistic MD simulations.

G.A.M.E.: GPU-accelerated mixture elucidator

GPU acceleration is useful in solving complex chemical information problems. Identifying unknown structures from the mass spectra of natural product mixtures has been a desirable yet unresolved issue in metabolomics. However, this elucidation process has been hampered by complex experimental data and the inability of instruments to completely separate different compounds. Fortunately, with current high-resolution mass spectrometry, one feasible strategy is to define this problem as extending a scaffold database with sidechains of different probabilities to match the high-resolution mass obtained from a high-resolution mass spectrum. By introducing a dynamic programming (DP) algorithm, it is possible to solve this NP-complete problem in pseudo-polynomial time. However, the running time of the DP algorithm grows by orders of magnitude as the number of mass decimal digits increases, thus limiting the boost in structural prediction capabilities. By harnessing the heavily parallel architecture of modern GPUs, we designed a “compute unified device architecture” (CUDA)-based GPU-accelerated mixture elucidator (G.A.M.E.) that considerably improves the performance of the DP, allowing up to five decimal digits for input mass data. As exemplified by four testing datasets with verified constitutions from natural products, G.A.M.E. allows for efficient and automatic structural elucidation of unknown mixtures for practical procedures.Graphical abstract

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.

Resolving the Conformational Dynamics of DNA with Ångstrom Resolution by Pulsed Electron-Electron Double Resonance and Molecular Dynamics

Pulsed electron-electron double resonance (PELDOR/DEER) experiments of nucleic acids with rigid spin labels provide highly accurate distance and orientation information. Here we combine PELDOR experiments with molecular dynamics (MD) simulations to arrive at an atomistic view of the conformational dynamics of DNA. The MD simulations closely reproduce the PELDOR time traces, and demonstrate that bending, in addition to twist-stretch motions, underpin the sub-μs dynamics of DNA. PELDOR experiments correctly rank DNA force fields and resolve subtle differences in the conformational ensembles of nucleic acids, on the order of 1-2 Å. Long-range distance and angle measurements with rigid spin labels provide critical input for the refinement of computer models and the elucidation of the structure and dynamics of complex biomolecules.

Understanding the mechanical response of double-stranded DNA and RNA under constant stretching forces using all-atom molecular dynamics

Multiple biological processes involve the stretching of nucleic acids (NAs). Stretching forces induce local changes in the molecule structure, inhibiting or promoting the binding of proteins, which ultimately affects their functionality. Understanding how a force induces changes in the structure of NAs at the atomic level is a challenge. Here, we use all-atom, microsecond-long molecular dynamics to simulate the structure of dsDNA and dsRNA subjected to stretching forces up to 20 pN. We determine all of the elastic constants of dsDNA and dsRNA and provide an explanation for three striking differences in the mechanical response of these two molecules: the threefold softer stretching constant obtained for dsRNA, the opposite twist-stretch coupling, and its nontrivial force dependence. The lower dsRNA stretching resistance is linked to its more open structure, whereas the opposite twist-stretch coupling of both molecules is due to the very different evolution of molecules' interstrand distance with the stretching force. A reduction of this distance leads to overwinding in dsDNA. In contrast, dsRNA is not able to reduce its interstrand distance and can only elongate by unwinding. Interstrand distance is directly correlated with the slide base-pair parameter and its different behavior in dsDNA and dsRNA traced down to changes in the sugar pucker angle of these NAs.

Probing the structure and in silico stability of cargo loaded DNA icosahedra using MD simulations

Platonic solids such as polyhedra based on DNA have been deployed for multifarious applications such as RNAi delivery, biological targeting and bioimaging. All of these applications hinge on the capability of DNA polyhedra for molecular display with high spatial precision. Therefore high resolution structural models of such polyhedra are critical to widen their applications in both materials and biology. Here, we present an atomistic model of a well-characterized DNA icosahedron, with demonstrated versatile functionalities in biological systems. We study the structure and dynamics of this DNA icosahedron using fully atomistic molecular dynamics (MD) simulation in explicit water and ions. The major modes of internal motion have been identified using principal component analysis. We provide a quantitative estimate of the radius of gyration (R_g), solvent accessible surface area (SASA) and volume of the icosahedron which is essential to estimate its maximal cargo carrying capacity. Importantly, our simulation of gold nanoparticles (AuNPs) encapsulated within DNA icosahedra revealed enhanced stability of the AuNP loaded DNA icosahedra compared to empty icosahedra. This is consistent with the experimental results that show high yields of cargo-encapsulated DNA icosahedra that have led to its diverse applications for precision targeting. These studies reveal that the stabilizing interactions between the cargo and the DNA scaffold powerfully position DNA polyhedra as targetable nanocapsules for payload delivery. These insights can be exploited for precise molecular display for diverse biological applications.

Long-range correlations in the mechanics of small DNA circles under topological stress revealed by multi-scale simulation

It is well established that gene regulation can be achieved through activator and repressor proteins that bind to DNA and switch particular genes on or off, and that complex metabolic networks determine the levels of transcription of a given gene at a given time. Using three complementary computational techniques to study the sequence-dependence of DNA denaturation within DNA minicircles, we have observed that whenever the ends of the DNA are constrained, information can be transferred over long distances directly by the transmission of mechanical stress through the DNA itself, without any requirement for external signalling factors. Our models combine atomistic molecular dynamics (MD) with coarse-grained simulations and statistical mechanical calculations to span three distinct spatial resolutions and timescale regimes. While they give a consensus view of the non-locality of sequence-dependent denaturation in highly bent and supercoiled DNA loops, each also reveals a unique aspect of long-range informational transfer that occurs as a result of restraining the DNA within the closed loop of the minicircles.

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.

Using Wavelet Analysis To Assist in Identification of Significant Events in Molecular Dynamics Simulations

Long time scale molecular dynamics (MD) simulations of biological systems are becoming increasingly commonplace due to the availability of both large-scale computational resources and significant advances in the underlying simulation methodologies. Therefore, it is useful to investigate and develop data mining and analysis techniques to quickly and efficiently extract the biologically relevant information from the incredible amount of generated data. Wavelet analysis (WA) is a technique that can quickly reveal significant motions during an MD simulation. Here, the application of WA on well-converged long time scale (tens of μs) simulations of a DNA helix is described. We show how WA combined with a simple clustering method can be used to identify both the physical and temporal locations of events with significant motion in MD trajectories. We also show that WA can not only distinguish and quantify the locations and time scales of significant motions, but by changing the maximum time scale of WA a more complete characterization of these motions can be obtained. This allows motions of different time scales to be identified or ignored as desired.

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.

Long-timescale dynamics of the Drew-Dickerson dodecamer

We present a systematic study of the long-timescale dynamics of the Drew–Dickerson dodecamer (DDD: d(CGCGAATTGCGC)₂) a prototypical B-DNA duplex. Using our newly parameterized PARMBSC1 force field, we describe the conformational landscape of DDD in a variety of ionic environments from minimal salt to 2 M Na⁺Cl⁻ or K⁺Cl⁻. The sensitivity of the simulations to the use of different solvent and ion models is analyzed in detail using multi-microsecond simulations. Finally, an extended (10 μs) simulation is used to characterize slow and infrequent conformational changes in DDD, leading to the identification of previously uncharacterized conformational states of this duplex which can explain biologically relevant conformational transitions. With a total of more than 43 μs of unrestrained molecular dynamics simulation, this study is the most extensive investigation of the dynamics of the most prototypical DNA duplex.