Molecular Dynamics Simulations of Nucleic Acids with a Generalized Born Solvation Model

Inclusion of Water Multipoles into the Implicit Solvation Framework Leads to Accuracy Gains

The current practical "workhorses" of the atomistic implicit solvation─the Poisson-Boltzmann (PB) and generalized Born (GB) models─face fundamental accuracy limitations. Here, we propose a computationally efficient implicit solvation framework, the Implicit Water Multipole GB (IWM-GB) model, that systematically incorporates the effects of multipole moments of water molecules in the first hydration shell of a solute, beyond the dipole water polarization already present at the PB/GB level. The framework explicitly accounts for coupling between polar and nonpolar contributions to the total solvation energy, which is missing from many implicit solvation models. An implementation of the framework, utilizing the GAFF force field and AM1-BCC atomic partial charges model, is parametrized and tested against the experimental hydration free energies of small molecules from the FreeSolv database. The resulting accuracy on the test set (RMSE ∼ 0.9 kcal/mol) is 12% better than that of the explicit solvation (TIP3P) treatment, which is orders of magnitude slower. We also find that the coupling between polar and nonpolar parts of the solvation free energy is essential to ensuring that several features of the IWM-GB model are physically meaningful, including the sign of the nonpolar contributions.

Improving the Accuracy of Physics-Based Hydration-Free Energy Predictions by Machine Learning the Remaining Error Relative to the Experiment

The accuracy of computational models of water is key to atomistic simulations of biomolecules. We propose a computationally efficient way to improve the accuracy of the prediction of hydration-free energies (HFEs) of small molecules: the remaining errors of the physics-based models relative to the experiment are predicted and mitigated by machine learning (ML) as a postprocessing step. Specifically, the trained graph convolutional neural network attempts to identify the "blind spots" in the physics-based model predictions, where the complex physics of aqueous solvation is poorly accounted for, and partially corrects for them. The strategy is explored for five classical solvent models representing various accuracy/speed trade-offs, from the fast analytical generalized Born (GB) to the popular TIP3P explicit solvent model; experimental HFEs of small neutral molecules from the FreeSolv set are used for the training and testing. For all of the models, the ML correction reduces the resulting root-mean-square error relative to the experiment for HFEs of small molecules, without significant overfitting and with negligible computational overhead. For example, on the test set, the relative accuracy improvement is 47% for the fast analytical GB, making it, after the ML correction, almost as accurate as uncorrected TIP3P. For the TIP3P model, the accuracy improvement is about 39%, bringing the ML-corrected model's accuracy below the 1 kcal/mol threshold. In general, the relative benefit of the ML corrections is smaller for more accurate physics-based models, reaching the lower limit of about 20% relative accuracy gain compared with that of the physics-based treatment alone. The proposed strategy of using ML to learn the remaining error of physics-based models offers a distinct advantage over training ML alone directly on reference HFEs: it preserves the correct overall trend, even well outside of the training set.

Distinct structure and gating mechanism in diverse NMDA receptors with GluN2C and GluN2D subunits

N-methyl-D-aspartate (NMDA) receptors are heterotetramers comprising two GluN1 and two alternate GluN2 (N2A-N2D) subunits. Here we report full-length cryo-EM structures of the human N1-N2D di-heterotetramer (di-receptor), rat N1-N2C di-receptor and N1-N2A-N2C tri-heterotetramer (tri-receptor) at a best resolution of 3.0 Å. The bilobate N-terminal domain (NTD) in N2D intrinsically adopts a closed conformation, leading to a compact NTD tetramer in the N1-N2D receptor. Additionally, crosslinking the ligand-binding domain (LBD) of two N1 protomers significantly elevated the channel open probability (Po) in N1-N2D di-receptors. Surprisingly, the N1-N2C di-receptor adopted both symmetric (minor) and asymmetric (major) conformations, the latter further locked by an allosteric potentiator, PYD-106, binding to a pocket between the NTD and LBD in only one N2C protomer. Finally, the N2A and N2C subunits in the N1-N2A-N2C tri-receptor display a conformation close to one protomer in the N1-N2A and N1-N2C di-receptors, respectively. These findings provide a comprehensive structural understanding of diverse function in major NMDA receptor subtypes.

The Intrinsic Radius as a Key Parameter in the Generalized Born Model to Adjust Protein-Protein Electrostatic Interaction

The generalized Born (GB) model is an extension of the continuum dielectric theory of Born solvation energy and is a powerful method for accelerating the molecular dynamic (MD) simulations of charged biological molecules in water. While the effective dielectric constant of water that varies as a function of the separation distance between solute molecules is incorporated into the GB model, adjustment of the parameters is indispensable for accurate calculation of the Coulomb (electrostatic) energy. One of the key parameters is the lower limit of the spatial integral of the energy density of the electric field around a charged atom, known as the intrinsic radius ρ. Although ad hoc adjustment of ρ has been conducted to improve the Coulombic (ionic) bond stability, the physical mechanism by which ρ affects the Coulomb energy remains unclear. Via energetic analysis of three differently sized systems, here, we clarify that the Coulomb bond stability increases with increasing ρ and that the increased stability is caused by the interaction energy term, not by the self-energy (desolvation energy) term, as was supposed previously. Our results suggest that the use of larger values for the intrinsic radii of hydrogen and oxygen atoms, together with the use of a relatively small value for the spatial integration cutoff in the GB model, can better reproduce the Coulombic attraction between protein molecules.

Molecular Insights into Substrate Binding of the Outer Membrane Enzyme OmpT

The enzyme OmpT of the outer membrane of Escherichia coli shows proteolytic activity and cleaves peptides and proteins. Using molecular dynamics simulations in a fully hydrated lipid bilayer on a time scale of hundreds of nanoseconds, we draw a detailed atomic picture of substrate recognition in the OmpT-holo enzyme complex. Hydrogen bonds and salt bridges are essential for maintaining the integrity of the active site and play a central role for OmpT in recognizing its substrate. Electrostatic interactions are critical at all stages from approaching the substrate to docking at the active site. Computational alanine scanning based on the Molecular Mechanics Generalized Born Surface Area (MM-GBSA) approach confirms the importance of multiple residues in the active site that form salt bridges. The substrate fluctuates along the axis of the β-barrel, which is associated with oscillations of the binding cleft formed by the residue pairs D210-H212 and D83-D85. Principal component analysis suggests that substrate and protein movements are correlated. We observe the transient presence of putative catalytic water molecules near the active site, which may be involved in the nucleophilic attack on the cleavable peptide bond of the substrate.

Molecular insight into the systematic affinity and selectivity of partner recognition sites between the WW1 and WW2 domains of human KIBRA neuroprotein

Human KIBRA, a member of the WWC family proteins, is an upstream regulator of the Salvador/Warts/Hippo (SWH) signaling pathway and predominately expressed in nervous system. The protein has two functionally regulatory domains WW1 and WW2 at N-terminal region, which recognize and bind to the PY-motif segments of their partner proteins to serve as a signaling scaffold role in the SWH pathway. The two domains are highly conserved, but their downstream ligands and biological functions may not be fully consistent. In this study, we attempted to systematically profile the PY-motif affinity to and selectivity between KIBRA WW1 and WW2 domains involved in partner recognition sites. Ontology mining was used to enrich the KIBRA-interacting proteins in literature libraries, from which a variety of PY-motif peptide segments were identified, and their binding behavior to each domain was then analyzed by integrating computational modeling and experimental assay. Most PY-motif peptides were found to interact potently with WW1 and WW2, but they generally only exhibit a moderate or modest selectivity between the two domains. Subsequently, several representative peptides were further examined in detail to elucidate the molecular mechanism underlying their affinity and selectivity. It is revealed that the middle motif region of PY-motif peptides is primarily responsible for the affinity and stability of peptide binding, but only contributes marginally to peptide selectivity. Instead, the N-terminal region and, particularly, C-terminal region of PY-motif peptides play a crucial role in the selectivity. Hydrophobic contacts and hydrogen bonds confer stability and specificity to the domain-peptide interaction, respectively.

A fast-slow method to treat solute dynamics in explicit solvent

Aiming to reduce the computational cost in the current explicit solvent molecular dynamics (MD) simulation, this paper proposes a fast-slow method for the fast MD simulation of biomolecules in explicit solvent. This fast-slow method divides the entire system into two parts: a core layer (typically solute or biomolecule) and a peripheral layer (typically solvent molecules). The core layer is treated using standard MD method but the peripheral layer is treated by a slower dynamics method to reduce the computational cost. We compared four different simulation models in testing calculations for several small proteins. These include gas-phase, implicit solvent, fast-slow explicit solvent and standard explicit solvent MD simulations. Our study shows that gas-phase and implicit solvent models do not provide a realistic solvent environment and fail to correctly produce reliable dynamic structures of proteins. On the other hand, the fast-slow method can essentially reproduce the same solvent effect as the standard explicit solvent model while gaining an order of magnitude in efficiency. This fast-slow method thus provides an efficient approach for accelerating the MD simulation of biomolecules in explicit solvent.

RNA-ligand molecular docking: advances and challenges

With rapid advances in computer algorithms and hardware, fast and accurate virtual screening has led to a drastic acceleration in selecting potent small molecules as drug candidates. Computational modeling of RNA-small molecule interactions has become an indispensable tool for RNA-targeted drug discovery. The current models for RNA-ligand binding have mainly focused on the docking-and-scoring method. Accurate docking and scoring should tackle four crucial problems: (1) conformational flexibility of ligand, (2) conformational flexibility of RNA, (3) efficient sampling of binding sites and binding poses, and (4) accurate scoring of different binding modes. Moreover, compared with the problem of protein-ligand docking, predicting ligand binding to RNA, a negatively charged polymer, is further complicated by additional effects such as metal ion effects. Thermodynamic models based on physics-based and knowledge-based scoring functions have shown highly encouraging success in predicting ligand binding poses and binding affinities. Recently, kinetic models for ligand binding have further suggested that including dissociation kinetics (residence time) in ligand docking would result in improved performance in estimating in vivo drug efficacy. More recently, the rise of deep-learning approaches has led to new tools for predicting RNA-small molecule binding. In this review, we present an overview of the recently developed computational methods for RNA-ligand docking and their advantages and disadvantages.

Calculation of Energy for RNA/RNA and DNA/RNA Duplex Formation by Molecular Dynamics Simulation

Abstract

The development of approaches for predictive calculation of hybridization properties of various nucleic acid (NA) derivatives is the basis for the rational design of the NA-based constructs. Modern advances in computer modeling methods provide the feasibility of these calculations. We have analyzed the possibility of calculating the energy of DNA/RNA and RNA/RNA duplex formation using representative sets of complexes (65 and 75 complexes, respectively). We used the classical molecular dynamics (MD) method, the MMPBSA or MMGBSA approaches to calculate the enthalpy (ΔH°) component, and the quasi-harmonic approximation (Q-Harm) or the normal mode analysis (NMA) methods to calculate the entropy (ΔS°) contribution to the Gibbs energy ($$\Delta G_{{37}}^{^\circ }$$ ) of the NA complex formation. We have found that the MMGBSA method in the analysis of the MD trajectory of only the NA duplex and the empirical linear approximation allow calculation of the enthalpy of formation of the DNA, RNA, and hybrid duplexes of various lengths and GC content with an accuracy of 8.6%. Within each type of complex, the combination of rather efficient MMGBSA and Q-Harm approaches being applied to the trajectory of only the bimolecular complex makes it possible to calculate the $$\Delta G_{{37}}^{^\circ }$$ of the duplex formation with an error value of 10%. The high accuracy of predictive calculation for different types of natural complexes (DNA/RNA, DNA/RNA, and RNA/RNA) indicates the possibility of extending the considered approach to analogs and derivatives of nucleic acids, which gives a fundamental opportunity in the future to perform rational design of new types of NA-targeted sequence-specific compounds.

Definition of the Binding Architecture to a Target Promoter of HP1043, the Essential Master Regulator of Helicobacter pylori

HP1043 is an essential orphan response regulator of Helicobacter pylori orchestrating multiple crucial cellular processes. Classified as a member of the OmpR/PhoB family of two-component systems, HP1043 exhibits a highly degenerate receiver domain and evolved to function independently of phosphorylation. Here, we investigated the HP1043 binding mode to a target sequence in the hp1227 promoter (P_hp1227). Scanning mutagenesis of HP1043 DNA-binding domain and consensus sequence led to the identification of residues relevant for the interaction of the protein with a target DNA. These determinants were used as restraints to guide a data-driven protein-DNA docking. Results suggested that, differently from most other response regulators of the same family, HP1043 binds in a head-to-head conformation to the P_hp1227 target promoter. HP1043 interacts with DNA largely through charged residues and contacts with both major and minor grooves of the DNA are required for a stable binding. Computational alanine scanning on molecular dynamics trajectory was performed to corroborate our findings. Additionally, in vitro transcription assays confirmed that HP1043 positively stimulates the activity of RNA polymerase.

Can Glycosylation Mask the Detection of MHC Expressing p53 Peptides by T Cell Receptors?

Proteins of the major histocompatibility complex (MHC) class I, or human leukocyte antigen (HLA) in humans interact with endogenous peptides and present them to T cell receptors (TCR), which in turn tune the immune system to recognize and discriminate between self and foreign (non-self) peptides. Of especial importance are peptides derived from tumor-associated antigens. T cells recognizing these peptides are found in cancer patients, but not in cancer-free individuals. What stimulates this recognition, which is vital for the success of checkpoint based therapy? A peptide derived from the protein p53 (residues 161–169 or p161) was reported to show this behavior. T cells recognizing this unmodified peptide could be further stimulated in vitro to create effective cancer killing CTLs (cytotoxic T lymphocytes). We hypothesize that the underlying difference may arise from post-translational glycosylation of p161 in normal individuals, likely masking it against recognition by TCR. Defects in glycosylation in cancer cells may allow the presentation of the native peptide. We investigate the structural consequences of such peptide glycosylation by investigating the associated structural dynamics.

Local Interactions That Contribute Minimal Frustration Determine Foldability

Earlier experiments suggest that the evolutionary information (conservation and coevolution) encoded in protein sequences is necessary and sufficient to specify the fold of a protein family. However, there is no computational work to quantify the effect of such evolutionary information on the folding process. Here we explore the role of early folding steps for sequences designed using coevolution and conservation through a combination of computational and experimental methods. We simulated a repertoire of native and designed WW domain sequences to analyze early local contact formation and found that the N-terminal β-hairpin turn would not form correctly due to strong non-native local contacts in unfoldable sequences. Through a maximum likelihood approach, we identified five local contacts that play a critical role in folding, suggesting that a small subset of amino acid pairs can be used to solve the "needle in the haystack" problem to design foldable sequences. Thus, using the contact probability of those five local contacts that form during the early stage of folding, we built a classification model that predicts the foldability of a WW sequence with 81% accuracy. This classification model was used to redesign WW domain sequences that could not fold due to frustration and make them foldable by introducing a few mutations that led to the stabilization of these critical local contacts. The experimental analysis shows that a redesigned sequence folds and binds to polyproline peptides with a similar affinity as those observed for native WW domains. Overall, our analysis shows that evolutionary-designed sequences should not only satisfy the folding stability but also ensure a minimally frustrated folding landscape.

Structure-Based Approaches to Classify the Functional Impact of ZBTB18 Missense Variants in Health and Disease

The Cys₂His₂ type zinc finger is a motif found in many eukaryotic transcription factor proteins that facilitates binding to genomic DNA so as to influence cellular gene expression. One such transcription factor is ZBTB18, characterized as a repressor that orchestrates the development of mammalian tissues including skeletal muscle and brain during embryogenesis. In humans, it has been recognized that disease-associated ZBTB18 missense variants mapping to the coding sequence of the zinc finger domain influence sequence-specific DNA binding, disrupt transcriptional regulation, and impair neural circuit formation in the brain. Furthermore, general population ZBTB18 missense variants that influence DNA binding and transcriptional regulation have also been documented within this domain; however, the molecular traits that explain why some variants cause disease while others do not are poorly understood. Here, we have applied five structure-based approaches to evaluate their ability to discriminate between disease-associated and general population ZBTB18 missense variants. We found that thermodynamic integration and Residue Scanning in the Schrodinger Biologics Suite were the best approaches for distinguishing disease-associated variants from general population variants. Our results demonstrate the effectiveness of structure-based approaches for the functional characterization of missense alleles to DNA binding, zinc finger transcription factor protein-coding genes that underlie human health and disease.

MD-TSPC4: Computational Method for Predicting the Thermal Stability of I-Motif

I-Motif is a tetrameric cytosine-rich DNA structure with hemi-protonated cytosine: cytosine base pairs. Recent evidence showed that i-motif structures in human cells play regulatory roles in the genome. Therefore, characterization of novel i-motifs and investigation of their functional implication are urgently needed for comprehensive understanding of their roles in gene regulation. However, considering the complications of experimental investigation of i-motifs and the large number of putative i-motifs in the genome, development of an in silico tool for the characterization of i-motifs in the high throughput scale is necessary. We developed a novel computation method, MD-TSPC4, to predict the thermal stability of i-motifs based on molecular modeling and molecular dynamic simulation. By assuming that the flexibility of loops in i-motifs correlated with thermal stability within certain temperature ranges, we evaluated the correlation between the root mean square deviations (RMSDs) of model structures and the thermal stability as the experimentally obtained melting temperature (Tm). Based on this correlation, we propose an equation for Tm prediction from RMSD. We expect this method can be useful for estimating the overall structure and stability of putative i-motifs in the genome, which can be a starting point of further structural and functional studies of i-motifs.

DNA-binding mechanisms of human and mouse cGAS: a comparative MD and MM/GBSA study

Cyclic GMP-AMP synthase (cGAS) can detect the presence of cytoplasmic DNA and activate the innate immune system via the cGAS-STING pathway. Although several structures of cGAS-DNA complexes were resolved recently, the molecular mechanism of cGAS in its recognition of DNA has not yet been fully understood. In order to reveal the subtle differences between human and mouse cGAS in terms of their DNA-binding mechanisms, four systems, both human and mouse cGAS in complex with two different DNA sequences of equal length, were studied by molecular dynamics simulations and molecular mechanics/generalized Born surface area analysis. Several residues, including ARG176/ARG161, ARG195/ARG180, ASN210/ASN196, LYS384/LYS372, CYM397/CYM385, LYS403/LYS391, LYS407/LYS395, and LYS411/LYS399, were identified to be the common key residues in the recognition of DNA for cGAS in both humans and mice. In addition, four residue pairs LYS173/ARG158, ASP177/LYS162, CYS199/LYS184, and GLU398/SER387 were suggested to be the major residues that make human cGAS and mouse cGAS different in terms of their binding to DNA. Besides the well-known zinc-thumb domain, two residues at the kink of the spine helix were also proposed for the first time to be the major binding motifs in cGAS-DNA interaction.

Multidimensional Global Optimization and Robustness Analysis in the Context of Protein-Ligand Binding

Accuracy of protein-ligand binding free energy calculations utilizing implicit solvent models is critically affected by parameters of the underlying dielectric boundary, specifically, the atomic and water probe radii. Here, a global multidimensional optimization pipeline is developed to find optimal atomic radii specifically for protein-ligand binding calculations in implicit solvent. The computational pipeline has these three key components: (1) a massively parallel implementation of a deterministic global optimization algorithm (VTDIRECT95), (2) an accurate yet reasonably fast generalized Born implicit solvent model (GBNSR6), and (3) a novel robustness metric that helps distinguish between nearly degenerate local minima via a postprocessing step of the optimization. A graph-based "kT-connectivity" approach to explore and visualize the multidimensional energy landscape is proposed: local minima that can be reached from the global minimum without exceeding a given energy threshold (kT) are considered to be connected. As an illustration of the capabilities of the optimization pipeline, we apply it to find a global optimum in the space of just five radii: four atomic (O, H, N, and C) radii and water probe radius. The optimized radii, ρW = 1.37 Å, ρC = 1.40 Å, ρH = 1.55 Å, ρN = 2.35 Å, and ρO = 1.28 Å, lead to a closer agreement of electrostatic binding free energies with the explicit solvent reference than two commonly used sets of radii previously optimized for small molecules. At the same time, the ability of the optimizer to find the global optimum reveals fundamental limits of the common two-dielectric implicit solvation model: the computed electrostatic binding free energies are still almost 4 kcal/mol away from the explicit solvent reference. The proposed computational approach opens the possibility to further improve the accuracy of practical computational protocols for binding free energy calculations.

Structural transitions in the RNA 7SK 5' hairpin and their effect on HEXIM binding

7SK RNA, as part of the 7SK ribonucleoprotein complex, is crucial to the regulation of transcription by RNA-polymerase II, via its interaction with the positive transcription elongation factor P-TEFb. The interaction is induced by binding of the protein HEXIM to the 5′ hairpin (HP1) of 7SK RNA. Four distinct structural models have been obtained experimentally for HP1. Here, we employ computational methods to investigate the relative stability of these structures, transitions between them, and the effects of mutations on the observed structural ensembles. We further analyse the results with respect to mutational binding assays, and hypothesize a mechanism for HEXIM binding. Our results indicate that the dominant structure in the wild type exhibits a triplet involving the unpaired nucleotide U40 and the base pair A43-U66 in the GAUC/GAUC repeat. This conformation leads to an open major groove with enough potential binding sites for peptide recognition. Sequence mutations of the RNA change the relative stability of the different structural ensembles. Binding affinity is consequently lost if these changes alter the dominant structure.

Truncated tetrahedral RNA nanostructures exhibit enhanced features for delivery of RNAi substrates

Using RNA as a material for nanoparticle construction provides control over particle size and shape at the nano-scale. RNA nano-architectures have shown promise as delivery vehicles for RNA interference (RNAi) substrates, allowing multiple functional entities to be combined on a single particle in a programmable fashion. Rather than employing a completely bottom-up approach to scaffold design, here multiple copies of an existing synthetic supramolecular RNA nano-architecture serve as building blocks along with additional motifs for the design of a novel truncated tetrahedral RNA scaffold, demonstrating that rationally designed RNA assemblies can themselves serve as modular pieces in the construction of larger rationally designed structures. The resulting tetrahedral scaffold displays enhanced characteristics for RNAi-substrate delivery in comparison to similar RNA-based scaffolds, as evidenced by its increased functional capacity, increased cellular uptake and ultimately an increased RNAi efficacy of its adorned Dicer substrate siRNAs. The unique truncated tetrahedral shape of the nanoparticle core appears to contribute to this particle's enhanced function, indicating the physical characteristics of RNA scaffolds merit significant consideration when designing platforms for delivery of functional RNAs via RNA nanoparticles.

Improved Poisson-Boltzmann Methods for High-Performance Computing

Implicit solvent models based on the Poisson-Boltzmann equation (PBE) have been widely used to study electrostatic interactions in biophysical processes. These models often treat the solvent and solute regions as high and low dielectric continua, leading to a large jump in dielectrics across the molecular surface which is difficult to handle. Higher order interface schemes are often needed to seek higher accuracy for PBE applications. However, these methods are usually very liberal in the use of grid points nearby the molecular surface, making them difficult to use on high-performance computing platforms. Alternatively, the harmonic average (HA) method has been used to approximate dielectric interface conditions near the molecular surface with surprisingly good convergence and is well suited for high-performance computing. By adopting a 7-point stencil, the HA method is advantageous in generating simple 7-banded coefficient matrices, which greatly facilitate linear system solution with dense data parallelism, on high-performance computing platforms such as a graphics processing unit (GPU). However, the HA method is limited due to its lower accuracy. Therefore, it would be of great interest for high-performance applications to develop more accurate methods while retaining the simplicity and effectiveness of the 7-point stencil discretization scheme. In this study, we have developed two new algorithms based on the spirit of the HA method by introducing more physical interface relations and imposing the discretized Poisson's equation to the second order, respectively. Our testing shows that, for typical biomolecules, the new methods significantly improve the numerical accuracy to that comparable to the second-order solvers and with ∼65% overall efficiency gain on widely available high-performance GPU platforms.

Thermodynamics of Helix formation in small peptides of varying lengthin vacuo, implicit solvent and explicit solvent: Comparison between AMBER force fields

Helix formation is of great significance in protein folding. The helix-forming tendencies of amino acids are accumulated along the sequence to determine the helix-forming tendency of peptides. Computer simulation can be used to model this process in atomic details and give structural insights. In the current work, we employ equilibrate-state free energy simulation to systematically study the folding/unfolding thermodynamics of a series of mutated peptides. Two AMBER force fields including AMBER99SB and AMBER14SB are compared. The new 14SB force field uses refitted torsion parameters compared with 99SB and they share the same atomic charge scheme. We find that in vacuo the helix formation is mutation dependent, which reflects the different helix propensities of different amino acids. In general, there are helix formers, helix indifferent groups and helix breakers. The helical structure becomes more favored when the length of the sequence becomes longer, which arises from the formation of additional backbone hydrogen bonds in the lengthened sequence. Therefore, the helix indifferent groups and helix breakers will become helix formers in long sequences. Also, protonation-dependent helix formation is observed for ionizable groups. In 14SB, the helical structures are more stable than in 99SB and differences can be observed in their grouping schemes, especially in the helix indifferent group. In solvents, all mutations are helix indifferent due to protein–solvent interactions. The decrease in the number of backbone hydrogen bonds is the same with the increase in the number of protein–water hydrogen bonds. The 14SB in explicit solvent is able to capture the free energy minima in the helical state while 14SB in implicit solvent, 99SB in explicit solvent and 99SB in implicit solvent cannot. The helix propensities calculated under 14SB agree with the corresponding experimental values, while the 99SB results obviously deviate from the references. Hence, implicit solvent models are unable to correctly describe the thermodynamics even for the simple helix formation in isolated peptides. Well-developed force fields and explicit solvents are needed to correctly describe the protein dynamics. Aside from the free energy, differences in conformational ensemble under different force fields in different solvent models are observed. The numbers of hydrogen bonds formed under different force fields agree and they are mostly determined by the solvent model.

Generalized Born Implicit Solvent Models for Biomolecules

It would often be useful in computer simulations to use an implicit description of solvation effects, instead of explicitly representing the individual solvent molecules. Continuum dielectric models often work well in describing the thermodynamic aspects of aqueous solvation and can be very efficient compared to the explicit treatment of the solvent. Here, we review a particular class of so-called fast implicit solvent models, generalized Born (GB) models, which are widely used for molecular dynamics (MD) simulations of proteins and nucleic acids. These approaches model hydration effects and provide solvent-dependent forces with efficiencies comparable to molecular-mechanics calculations on the solute alone; as such, they can be incorporated into MD or other conformational searching strategies in a straightforward manner. The foundations of the GB model are reviewed, followed by examples of newer, emerging models and examples of important applications. We discuss their strengths and weaknesses, both for fidelity to the underlying continuum model and for the ability to replace explicit consideration of solvent molecules in macromolecular simulations.

Identification of AtHsp90.6 involved in early embryogenesis and its structure prediction by molecular dynamics simulations

Heat-shock protein of 90 kDa (Hsp90) is a key molecular chaperone involved in folding the synthesized protein and controlling protein quality. Conformational dynamics coupled to ATPase activity in N-terminal domain is essential for Hsp90's function. However, the relevant process is still largely unknown in plant Hsp90s, especially those required for plant embryogenesis which is inextricably tied up with human survival. Here, AtHsp90.6, a member of Hsp90 family in Arabidopsis, was firstly identified as a protein essential for embryogenesis. Thus we modelled AtHsp90.6 in its functionally closed 'lid-down' and open 'lid-up' states, exploring the nucleotide binding mechanism in these two states. Free energy landscape and electrostatic potential analysis revealed the switching mechanism between these two states. Collectively, this study quantitatively analysed the conformational changes of AtHsp90.6 bound to ATP or ADP. This result may help us understand the mechanism of action of AtHsp90.6 in future.

Systematic analysis and identification of unexpected interactions from the neuroprotein drug interactome in hydrocephalus pharmacological intervention

Hydrocephalus is a neurological condition caused by an abnormal accumulation of cerebrospinal fluid; pharmacological intervention of the disease has been found to elicit a variety of adverse drug reactions (ADRs) in central nervous system (CNS) by unexpectedly targeting certain functional neuroproteins. Here, a systematic neuroprotein drug interactome (SNDI) is created for 11 hydrocephalus drugs/metabolites plus 20 control drugs across 518 druggable pockets on the surface of 472 CNS neuroproteins via a large-scale molecular docking approach. Heuristic clustering analysis of the SNDI profile divides the 31 investigated drug ligands into a distinct panel and a background panel; the former consists of two hydrocephalus drugs (Furosemide and Triamterene) and their respective metabolites (Furosemide glucuronide and Hydroxytriamterene) that are inferred to have generally high affinity towards the whole array of neuroprotein pockets. A total of 13 neuroproteins are enriched in gene ontology semantic mining as putative unexpected targets of the distinct panel, and their intermolecular interactions with hydrocephalus drugs/metabolites are investigated in detail using dynamics simulation and energetics analysis. We also perform kinase assay and viability test to substantiate the interactome analysis. It is found that the Furosemide and Triamterene have significant cytotoxic effects on normal human astrocytes, in which the Triamterene can inhibit the neurokinase ROCK2, a representative of putative unexpected targets, with a high activity, which is comparable with the sophisticated ROCK2 inhibitor Fasudil.

An efficient second-order poisson-boltzmann method

Immersed interface method (IIM) is a promising high-accuracy numerical scheme for the Poisson-Boltzmann model that has been widely used to study electrostatic interactions in biomolecules. However, the IIM suffers from instability and slow convergence for typical applications. In this study, we introduced both analytical interface and surface regulation into IIM to address these issues. The analytical interface setup leads to better accuracy and its convergence closely follows a quadratic manner as predicted by theory. The surface regulation further speeds up the convergence for nontrivial biomolecules. In addition, uncertainties of the numerical energies for tested systems are also reduced by about half. More interestingly, the analytical setup significantly improves the linear solver efficiency and stability by generating more precise and better-conditioned linear systems. Finally, we implemented the bottleneck linear system solver on GPUs to further improve the efficiency of the method, so it can be widely used for practical biomolecular applications. © 2019 Wiley Periodicals, Inc.

Structural Analysis of a Trimer of β2-Microgloblin Fragment by Molecular Dynamics Simulations

A peptide β₂-m_21-31, which is a fragment from residue 21 to residue 31 of β₂-microgloblin, is experimentally known to self-assemble and form amyloid fibrils. In order to understand the mechanism of amyloid fibril formations, we applied the replica-exchange molecular dynamics method to the system consisting of three fragments of β₂-m_21-31. From the analyses on the temperature dependence, we found that there is a clear phase transition temperature in which the peptides aggregate with each other. Moreover, we found by the free energy analyses that there are two major stable states: One of them is like amyloid fibrils and the other is amorphous aggregates.

Combining Structural Data with Computational Methodologies to Investigate Structure-Function Relationships in TRP Channels

Since the emergence of high-resolution three-dimensional structures of membrane proteins, and the increasing availability of state-of-the-art algorithms and high-performance-computing facilities, classical molecular dynamics (MD) simulations have become a routine device to explore the molecular behavior of these proteins. The rise of cryo-electron microscopy (cryo-EM) as a credible experimental tool to resolve structures at an atomic level has revolutionized structural biology in recent years, culminating in the disclosure of the first high-resolution three-dimensional structure of a transient receptor potential (TRP) channel, the vanilloid receptor 1 (TRPV1). As a result, the number of research articles investigating the molecular behavior of TRP channels using macromolecular simulation techniques has proliferated. This review provides an overview of the current state of this field, including our understanding of TRP channel structure, the framework of classical MD simulations, and how to perform such simulations to investigate structure-function relationships in TRP channels.

Two-Photon Absorption Properties of Gold Fluorescent Protein: A Combined Molecular Dynamics and Quantum Chemistry Study

Molecular dynamics (MD) simulations were carried out to obtain the conformational changes of the chromophore in the gold fluorescent protein (PDB ID: 1OXF ). To obtain two-photon absorption (TPA) cross-sections, time dependent density functional theory (TD-DFT) computations were performed for chromophore geometries sampled along the trajectory. The TD-DFT computations used the CAM-B3LYP functional and 6-31+G(d) basis set. Results showed that two dihedral angles change remarkably over the simulation time. TPA cross-sections were found to average 13.82 GM for the excitation to S₁ computed from the equilibrium geometries; however, extending the structures with a water molecule and GLU residue, which make H bonds with the chromophore molecule, increased excitation energies and TPA cross-sections significantly. Besides the effects of the surrounding residues and the dihedrals on the spectroscopic properties, some bond lengths affected the excitation energies and the TPA cross-sections significantly (up to ±25-30%), while the effects of the bond angles were smaller (±5%). Overall the present results provide insight into the effects of the conformational flexibility on TPA (with gold fluorescent protein as a specific example) and suggest that further experimental measurements of TPA for the gold fluorescent protein should be undertaken.

In silico-in vitro discovery of untargeted kinase-inhibitor interactions from kinase-targeted therapies: A case study on the cancer MAPK signaling pathway

Protein kinase inhibitors have been widely used as therapeutic agents to treat a variety of diseases, but many of them may cause off-target effects by unexpectedly targeting other noncognate kinases due to high conversion across the protein kinase family. The mitogen-activated protein kinase (MAPK) signaling pathway plays an essential role in tumorigenesis, which has been recognized as a high priority in the druggable target candidates of anticancer therapy. Here, we attempt to investigate the untargeted kinase-inhibitor interactions (UKIIs) of kinase-targeted therapies for the cancer MAPK signaling cascade via an integration of biomolecular modeling, cell viability assay and kinase inhibition analysis. A systematic kinase-inhibitor interaction profile is created for 28 FDA-approved kinase inhibitor drugs across 9 caner-related MAPK kinases. The created profile is analyzed at structural, energetic and dynamic levels and, consequently, totally 18 promising UKII pairs with high theoretical affinity are derived, from which the noncognate inhibitors Cabozantinib, Regorafenib and Crizotinib are selected to test their cytotoxic effects on human epithelial colorectal adenocarcinoma Caco-2 cell line and inhibition activity against the recombinant protein of human p38α kinase domain. The obtained results are compared with two cognate MAPK inhibitors JNK-IN-8 and BIRB796. As might be expected, the Regorafenib, Crizotinib and Cabozantinib exhibit high, moderate and low cytotoxicities, respectively. In addition, the Regorafenib is determined to have a potent p38α-inhibitory activity. This is basically in line with the test results of positive controls JNK-IN-8 and BIRB796 and can be well confirmed by computational modeling.

Interaction of microtubule depolymerizing agent indanocine with different human αβ tubulin isotypes

Tubulin isotypes are known to regulate the stability and dynamics of microtubules, and are also involved in the development of resistance against microtubule-targeted cancer drugs. Indanocine, a potent microtubule depolymerizing agent, is highly active against multidrug-resistant (MDR) cancer cells without affecting normal cells. It is known to disrupt microtubule dynamics in cells and induce apoptotic cell death. Indanocine is reported to bind to tubulin at the colchicine site i.e. at the interface of αβ tubulin heterodimer. However, it’s precise binding mode, involved molecular interactions and the binding affinities with different αβ-tubulin isotypes present in MDR cells are not well understood. Here, the binding affinities of human αβ-tubulin isotypes with indanocine were examined, employing the molecular modeling approach i.e. docking, molecular dynamics simulation and binding energy calculations. Multiple sequence analysis suggests that the amino acid sequences are different in the indanocine binding pockets of βI, βIIa, βIII and βVI isotypes. However, such differences are not observed in the amino acid sequences of βIVa, βIVb, and βV tubulin isotypes at indanocine binding pockets. Docking and molecular dynamics simulation results show that indanocine prefers the interface binding pocket of αβIIa, αβIII, αβIVb, αβV, and αβVI tubulin isotypes; whereas it is expelled from the interface binding pocket of αβIVa and αβI-tubulin isotypes. Further, binding free energy calculations show that αβVI has the highest binding affinity and αβI has the lowest binding affinity for indanocine among all β-tubulin isotypes. The binding free energy decreases in the order of αβVI > αβIVb > αβIIa > αβIII > αβV > αβIVa > αβI. Thus, our study provides a significant understanding of involved molecular interactions of indanocine with tubulin isotypes, which may help to design potent indanocine analogues for specific tubulin isotypes in MDR cells in future.

A molecular mechanics study on GA codon box translation

The GA codon box incorporates the two-fold degeneracy of aspartic acid and of glutamic acid. Using the molecular mechanics approach of the AMBER suite, the four codons of the GA box are paired via H-bonding with two aspartic acid anticodons and two glutamic acid anticodons to yield 8 cognate and 11 non-cognate codon-anticodon duplexes. In addition four select non-cognate duplexes between the GA box codons and three alanine anticodons are also studied. These 23 duplexes display a variety of base-pairing possibilities at the wobble position. Cognate duplexes are differentiated from non-cognate duplexes on the grounds of structure and stability (chiefly the former). The results are in line with Crick's wobble hypothesis, and corroborate the observed reading properties of the aspartic acid anticodons GUC and QUC and of the glutamic acid anticodons CUC and S^mnUC.

Identification and characterization of the druggable kinase targets of olmesartan and its analogues from a systematic kinase-chemical interaction profile in atherosclerosis

Olmesartan (OL) is the pharmacologically active metabolite of Olmesartan medoxomil (OM), an FDA-approved angiotensin II receptor antagonist for administrating cardiovascular diseases. The drug has been found to have potential effects on diverse protein kinase signaling involved in the pathogenesis of atherosclerosis, either by directly inhibiting the hub kinases or by indirectly modulating marginal members in the signaling pathways. In the present study, we computationally model the kinase-chemical Interaction Profile between six OL-related chemicals (i.e. OL, OM, Valsartan [VL], Losartan [LS], Candesartan [CD] and Telmisartan [TL]) and 23 human protein kinases in atherosclerosis. The profile is analyzed systematically at molecular level to identify unexpected kinase targets for OL. There is a good consistence between co-citation frequency and affinity scoring for the chemical association with kinase candidates; the OL and its analogs VL and LS exhibit a similar binding profile to the atherosclerosis kinase spectrum. It is suggested that the Ser/Thr-specific kinases PI3Kα and ROCK1 are potential druggable targets of OL for atherosclerosis therapy. As a paradigm, kinase assays reveal that the inhibitory potency of OL and Y-27632 (positive control) on ROCK1 is determined at micromolar level, while the OM (negative control) possesses no detectable activity for the kinase.

Accuracy Comparison of Generalized Born Models in the Calculation of Electrostatic Binding Free Energies

The need for accurate yet efficient representation of the aqueous environment in biomolecular modeling has led to the development of a variety of generalized Born (GB) implicit solvent models. While many studies have focused on the accuracy of available GB models in predicting solvation free energies, a systematic assessment of the quality of these models in binding free energy calculations, crucial for rational drug design, has not been undertaken. Here, we evaluate the accuracies of eight common GB flavors (GB-HCT, GB-OBC, GB-neck2, GBNSR6, GBSW, GBMV1, GBMV2, and GBMV3), available in major molecular dynamics packages, in predicting the electrostatic binding free energies ( ΔΔ G_el) for a diverse set of 60 biomolecular complexes belonging to four main classes: protein-protein, protein-drug, RNA-peptide, and small complexes. The GB flavors are examined in terms of their ability to reproduce the results from the Poisson-Boltzmann (PB) model, commonly used as accuracy reference in this context. We show that the agreement with the PB of ΔΔ G_el estimates varies widely between different GB models and also across different types of biomolecular complexes, with R² correlations ranging from 0.3772 to 0.9986. A surface-based "R6" GB model recently implemented in AMBER shows the closest overall agreement with reference PB ( R² = 0.9949, RMSD = 8.75 kcal/mol). The RNA-peptide and protein-drug complex sets appear to be most challenging for all but one model, as indicated by the large deviations from the PB in ΔΔ G_el. Small neutral complexes present the least challenge for most of the GB models tested. The quantitative demonstration of the strengths and weaknesses of the GB models across the diverse complex types provided here can be used as a guide for practical computations and future development efforts.

Electrostatic Interactions in Protein Structure, Folding, Binding, and Condensation

Charged and polar groups, through forming ion pairs, hydrogen bonds, and other less specific electrostatic interactions, impart important properties to proteins. Modulation of the charges on the amino acids, e.g., by pH and by phosphorylation and dephosphorylation, have significant effects such as protein denaturation and switch-like response of signal transduction networks. This review aims to present a unifying theme among the various effects of protein charges and polar groups. Simple models will be used to illustrate basic ideas about electrostatic interactions in proteins, and these ideas in turn will be used to elucidate the roles of electrostatic interactions in protein structure, folding, binding, condensation, and related biological functions. In particular, we will examine how charged side chains are spatially distributed in various types of proteins and how electrostatic interactions affect thermodynamic and kinetic properties of proteins. Our hope is to capture both important historical developments and recent experimental and theoretical advances in quantifying electrostatic contributions of proteins.

Proteomic analysis of stachyose contribution to the growth of Lactobacillus acidophilus CICC22162

Stachyose is a functional oligosaccharide, acting as a potential prebiotic for colonic fermentation.

Solid-State NMR Studies of Amyloid Materials: A Protocol to Define an Atomic Model of Aβ(1-42) in Amyloid Fibrils

Intense efforts have been made to understand the molecular structures of misfolded amyloid β (Aβ) in order to gain insight into the pathological mechanism of Alzheimer's disease. Solid-state NMR spectroscopy (SSNMR) is considered a primary tool for elucidating the structures of insoluble and noncrystalline amyloid fibrils and other amyloid assemblies. In this chapter, we describe a detailed protocol to obtain the first atomic model of the 42-residue human Aβ peptide Aβ(1-42) in structurally homogeneous amyloid fibrils from our recent SSNMR study (Nat Struct Mol Biol 22:499-505, 2015). Despite great biological and clinical interest in Aβ(1-42) fibrils, their structural details have been long-elusive until this study. The protocol is divided into four sections. First, the solid-phase peptide synthesis (SPPS) and purification of monomeric Aβ(1-42) is described. We illustrate a controlled incubation method to prompt misfolding of Aβ(1-42) into homogeneous amyloid fibrils in an aqueous solution with fragmented Aβ(1-42) fibrils as seeds. Next, we detail analysis of Aβ(1-42) fibrils by SSNMR to obtain structural restraints. Finally, we describe methods to construct atomic models of Aβ(1-42) fibrils based on SSNMR results through two-stage molecular dynamics calculations.

Molecular Modeling of Multidrug Properties of Resistance Nodulation Division (RND) Transporters

Efflux pumps of the resistance nodulation division (RND) superfamily are among the major contributors to intrinsic and acquired multidrug resistance in Gram-negative bacteria. Structural information on AcrAB-TolC and MexAB-OprM, major efflux pumps of Escherichia coli and Pseudomonas aeruginosa respectively, boosted intensive research aimed at understanding the molecular mechanisms ruling the active extrusion processes. In particular, several studies were devoted to the understanding of the determinants behind the extraordinary broad specificity of the RND transporters AcrB and MexB. In this chapter, we discuss the ever-growing role computational methods have been playing in deciphering key structural and dynamical features of these transporters and of their interaction with substrates and inhibitors. We further discuss and illustrate examples from our lab of how molecular docking, homology modeling, all-atom molecular dynamics simulations and in silico free energy estimations can all together give precious insights into the processes of recognition and extrusion of substrates, as well as on the possible inhibition strategies.

Molecular insight into mutation-induced conformational change in metastasic bowel cancer BRAF kinase domain and its implications for selective inhibitor design

Oncogenic BRAF V600E mutation confers constitutive activation for the kinase and is closely related to the pathogenesis of metastasic bowel cancer (MBC). Here, the V600E-induced conformational change in MBC BRAF kinase domain is characterized systematically at structural, energetic and dynamic levels. The mutation is observed to cause a conformational conversion of the kinase's activation loop from DFG-out to DFG-in, thus activating the kinase. Electrostatic force is primarily responsible for the conformational conversion and stabilization of DFG-in associated with the mutation. Molecular docking calculations are employed to analyze the binding mode difference of mutant-selective inhibitors between the DFG-out and DFG-in conformations of BRAF kinase. It is revealed that the mutation can reshape inhibitor selectivity profile by altering kinase loop conformation. Several compounds are determined to have a high or moderate selectivity for mutant over wild-type kinase. The selectivity is primarily originated from hydrogen bond interactions of inhibitor ligands with mutant rather than wild type due to the conformational difference in kinase domain.

Conformational rearrangement of 1,2-d(GG) intrastrand cis-diammineplatinum crosslinked DNA is driven by counter-ion penetration within the minor groove of the modified site

The major structural aberrations of DNA induced by a cis-diammineplatinum (II) 1,2-d(GG) intrastrand cross-link (CPT) have been known for decades. To gain deeper insights into the structural dynamics of the sequence-dependent DNA distortions adjacent to the CPT adduct, we employed molecular modeling and molecular dynamics (MD) simulations. The structural dynamics of native (N-DNA) and cisPt 1,2-d(GG) crosslinked (CPT-DNA) in the form of symmetric 36 nt d(G₂T₁₅G*G*T₁₅G₂)●C₂A₁₅CCA₁₅C₂) oligonucleotide duplexes is compared. The selected sequence context enabled tracking of the origin of the DNA axis curvature at the YpR flexible points (N-DNA), the enhancement of axis bending, and further distortions due to steric/electrostatic perturbations arising from the CPT-crosslink. In addition to the known structural distortions of CPT-DNA: helix bend towards the major groove; local helix unwinding; high roll angle between cross-linked guanine bases; and adoption of A-form DNA on the 5'-side of the CPT-crosslink (TpG junction); our results show the existence of a singular irreversible and reproducible conformational rearrangement, not previously observed, resulting in two stable CPT-DNA1 and CPT-DNA2 conformers. The CPT-DNA2 conformation presents an enhanced DNA axis bend and a wider and shallower minor grove with increased solvent accessibility within the modified site. It is concluded that the polymorphous (unstable) DNA environment near the cisPt 1,2-d(GG) unit in synergy with specific dynamic events, such as prolonged minor groove retention of particular Na⁺ ions and water redistribution within the d(TG*G*T) site, together with the formation of extra and more stable H-bonds between Pt(NH₃)₂ amines and neighboring nucleotides, are cooperatively responsible for the initiation of the conformational rearrangement leading to the CPT-DNA2 conformer, which, surprisingly, closely resembles the HMGB1-bound CPT-DNA structure. Graphical abstract Superimposed averaged structures of normal (N-DNA, green) and cisplatin intrastrand cross-linked (CPT-DNA, orange). Global DNA axes: N-DNA (blue beads); CPT-DNA (red beads); PT (yellow sphere); crosslinked dGs viewed from the minor groove (bold).

Integrated discovery of FOXO1-DNA stabilizers from marine natural products to restore chemosensitivity to anti-EGFR-based therapy for metastatic lung cancer

The transcription factor forkhead box O1 (FOXO1) negatively regulates activated EGFR signaling by turning on the gene expression of tumor suppressor Kruppel-like factor 6. Here, we propose that the chemosensitivity to anti-EGFR-based lung cancer therapy can be restored by stabilization of the FOXO1-DNA complex architecture using small-molecule marine natural medicines. A synthetic protocol that integrates computational ligand-protein-DNA binding analysis and an experimental fluorescence binding assay was applied against a large library of structurally diverse, drug-like marine natural products to discover novel stabilizers of DNA-bound FOXO1 conformation. The screening utilized chemical similarity analysis to exclude structurally redundant compounds, and then carried out high-throughput molecular docking and computational binding analysis to identify potential marine natural product candidates. Consequently, eight commercially available hits were selected and tested in vitro, from which four marine natural product compounds (tanzawaic acid D, hymenidin, cribrostatin 6 and barbamide) were found to have high or moderate potency to selectively bind to the FOXO1 DNA-binding domain (DBD) in the presence of its cognate DNA partner. Atomistic molecular dynamics (MD) simulations revealed that the identified stabilizers do not directly interact with DNA; instead, they can effectively stabilize the free FOXO1 DBD domain in the DNA-bound conformation and thus promote the binding of FOXO1 to DNA.

Over-Destabilization of Protein-Protein Interaction in Generalized Born Model and Utility of Energy Density Integration Cutoff

The generalize Born (GB) model is frequently used in MD simulations of biomolecular systems in aqueous solution. The GB model is usually based on the so-called Coulomb field approximation (CFA) for the energy density integration. In this study, we report that the GB model with CFA overdestabilizes the long-range electrostatic attraction between oppositely charged molecules (ionic bond forming two-helix system and kinesin-tubulin system) when the energy density integration cutoff, r_max, which is used to calculate the Born energy, is set to a large value. We show that employing large r_max, which is usually expected to make simulation results more accurate, worsens the accuracy so that the attraction is changed into repulsion. It is demonstrated that the overdestabilization is caused by the overestimation of the desolvation penalty upon binding that originates from CFA. We point out that the overdestabilization can be corrected by employing a relatively small cutoff (r_max = 10-15 Å), affirming that the GB models, even with CFA, can be used as a powerful tool to theoretically study the protein-protein interaction, particularly on its dynamical aspect, such as binding and unbinding.

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.

Folding behavior of a T-shaped, ribosome-binding translation enhancer implicated in a wide-spread conformational switch

Turnip crinkle virus contains a T-shaped, ribosome-binding, translation enhancer (TSS) in its 3'UTR that serves as a hub for interactions throughout the region. The viral RNA-dependent RNA polymerase (RdRp) causes the TSS/surrounding region to undergo a conformational shift postulated to inhibit translation. Using optical tweezers (OT) and steered molecular dynamic simulations (SMD), we found that the unusual stability of pseudoknotted element H4a/Ψ3 required five upstream adenylates, and H4a/Ψ3 was necessary for cooperative association of two other hairpins (H5/H4b) in Mg2+. SMD recapitulated the TSS unfolding order in the absence of Mg2+, showed dependence of the resistance to pulling on the 3D orientation and gave structural insights into the measured contour lengths of the TSS structure elements. Adenylate mutations eliminated one-site RdRp binding to the 3'UTR, suggesting that RdRp binding to the adenylates disrupts H4a/Ψ3, leading to loss of H5/H4b interaction and promoting a conformational switch interrupting translation and promoting replication.

Direct binding of microRNA-21 pre-element with Regorafenib: An alternative mechanism for anti-colorectal cancer chemotherapy?

The Regorafenib is a broad-spectrum kinase inhibitor that has been approved to treat colorectal cancer (CRC). However, evidences have shown that the agent is also implicated in drug interaction with microRNA-21 (miR-21), an oncogenic miRNA which plays a key role in resisting programmed cell death in CRC cells. Here, we supposed that, instead of kinase inhibition, Regorafenib can directly bind to and then stabilize miR-21 pre-element, thus preventing RNase Dicer-meditated cleavage of the pre-element to mature miR-21. In order to verify the notion, an in silico-in vitro integrated investigation of the direct intermolecular interaction between Regorafenib and miR-21 pre-element was performed by using active pocket identification, RNA-ligand docking, molecular dynamics (MD) simulation, binding energetic analysis, and fluorescence-based assay. It was revealed that the Regorafenib can bind at the major groove-like stem region of miR-21 pre-element through three geometrically satisfactory hydrogen bonds (H-bonds) as well as a number of hydrophobic forces and π-π stacking, conferring strong specificity and high stability to the RNA-ligand complex system (K_d=0.73μM). Separate inversion mutation of two base pairs (G6C, C12G) and (A13U, U4A) that are involved in the H-bonding can considerably impair the affinity of Regorafenib to miR-21 pre-element, with K_d increase to 27 and 96μM, respectively. All these supported that Regorafenib can directly bind to miR-21 pre-element at molecular level and the binding mode can be properly modeled by using the proposed integrated strategy. This study would provide a potential, alternative mechanism for anti-colorectal cancer chemotherapy with Regorafenib.

Histone Acetylation Regulates Chromatin Accessibility: Role of H4K16 in Inter-nucleosome Interaction

The N-terminal tail of histone H4 is an indispensable mediator for inter-nucleosome interaction, which is required for chromatin fiber condensation. H4K16 acetylation (H4K16Ac) activates gene transcription by influencing both chromatin structure and interplay with nonhistone proteins. To understand the influence of H4K16Ac on inter-nucleosome interaction, we performed a simulation study for the H4 tail in the context of two nucleosomes in neighboring unit cells in the crystal structure. The binding conformation of H4 tail with/without K16Ac was sampled by replica exchange with solute tempering, and the free energy landscape was explored by metadynamics. The results indicate two important features of H4K16: 1) it is the first button to anchor the H4 tail on the adjacent nucleosome; and 2) it is the only acetylation site interacting with the acidic patch. H4K16Ac disrupts the electrostatic interactions of K16, weakens H4 tail-acidic patch binding, and significantly increases H4 tail conformation diversity. Our study suggests that H4K16Ac directly reduces the inter-nucleosome interaction mediated by the H4 tail, which might further encourage the binding of nonhistone proteins on the acidic patch.