Evaluation of the stereochemical quality of predicted RNA 3D models in the RNA-Puzzles submissions.

In silico prediction is a well-established approach to derive a general shape of an RNA molecule based on its sequence or secondary structure. This paper reports an analysis of the stereochemical quality of the RNA three-dimensional models predicted using dedicated computer programs. The stereochemistry of 1052 RNA 3D structures, including 1030 models predicted by fully automated and human-guided approaches within 22 RNA-Puzzles challenges and reference structures, is analyzed. The evaluation is based on standards of RNA stereochemistry that the Protein Data Bank requires from deposited experimental structures. Deviations from standard bond lengths and angles, planarity, or chirality are quantified. A reduction in the number of such deviations should help in the improvement of RNA 3D structure modeling approaches.

Advanced Sampling Methods for Multiscale Simulation of Disordered Proteins and Dynamic Interactions.

Intrinsically disordered proteins (IDPs) are highly prevalent and play important roles in biology and human diseases. It is now also recognized that many IDPs remain dynamic even in specific complexes and functional assemblies. Computer simulations are essential for deriving a molecular description of the disordered protein ensembles and dynamic interactions for a mechanistic understanding of IDPs in biology, diseases, and therapeutics. Here, we provide an in-depth review of recent advances in the multi-scale simulation of disordered protein states, with a particular emphasis on the development and application of advanced sampling techniques for studying IDPs. These techniques are critical for adequate sampling of the manifold functionally relevant conformational spaces of IDPs. Together with dramatically improved protein force fields, these advanced simulation approaches have achieved substantial success and demonstrated significant promise towards the quantitative and predictive modeling of IDPs and their dynamic interactions. We will also discuss important challenges remaining in the atomistic simulation of larger systems and how various coarse-grained approaches may help to bridge the remaining gaps in the accessible time- and length-scales of IDP simulations.

Computational and biophysical determination of pathogenicity of variants of unknown significance in cardiac thin filament.

Point mutations within sarcomeric proteins have been associated with altered function and cardiomyopathy development. Difficulties remain, however, in establishing the pathogenic potential of individual mutations, often limiting the use of genotype in management of affected families. To directly address this challenge, we utilized our all-atom computational model of the human full cardiac thin filament (CTF) to predict how sequence substitutions in CTF proteins might affect structure and dynamics on an atomistic level. Utilizing molecular dynamics calculations, we simulated 21 well-defined genetic pathogenic cardiac troponin T and tropomyosin variants to establish a baseline of pathogenic changes induced in computational observables. Computational results were verified via differential scanning calorimetry on a subset of variants to develop an experimental correlation. Calculations were performed on 9 independent variants of unknown significance (VUS), and results were compared with pathogenic variants to identify high-resolution pathogenic signatures. Results for VUS were compared with the baseline set to determine induced structural and dynamic changes, and potential variant reclassifications were proposed. This unbiased, high-resolution computational methodology can provide unique structural and dynamic information that can be incorporated into existing analyses to facilitate classification both for de novo variants and those where established approaches have provided conflicting information.

Updating atomic charge parameters of aliphatic amino acids: a quest to improve the performance of molecular modeling via sequential molecular dynamics and DFT-GIAO-NMR calculations.

In this work, we observe the behavior of the dipole moment, atomic charges, solute-solvent interactions and NMR spectroscopy of aliphatic amino acids in a water solution via the computational simulations of classical molecular dynamics and DFT quantum calculations. Our results indicate that the convergence of the atomic charge of the solute, from an iterative process, together with the dipole moment of the amino acid, alters the lifetime of hydrogen bonds present in the first solvation shell, resulting in the modification of its structure and dynamics. Using GIAO-DFT-NMR calculations, we assessed the impact of these structural solute-solvent modifications on the magnetic shielding constants of the solute carbon atoms. In this sense, we evaluate the importance of an update in parameters that describe atomic charges present in the CHARMM36 force field.

Estimating computational limits on theoretical descriptions of biological cells.

There has been much success recently in theoretically simulating parts of complex biological systems on the molecular level, with the goal of first-principles modeling of whole cells. However, there is the question of whether such simulations can be performed because of the enormous complexity of cells. We establish approximate equations to estimate computation times required to simulate highly simplified models of cells by either molecular dynamics calculations or by solving molecular kinetic equations. Our equations place limits on the complexity of cells that can be theoretically understood with these two methods and provide a first step in developing what can be considered biological uncertainty relations for molecular models of cells. While a molecular kinetics description of the genetically simplest bacterial cell may indeed soon be possible, neither theoretical description for a multicellular system, such as the human brain, will be possible for many decades and may never be possible even with quantum computing.

Design, synthesis, and evaluation of "dual-site"-binding diarylpyrimidines targeting both NNIBP and the NNRTI adjacent site of the HIV-1 reverse transcriptase.

Inspired by our previous efforts to improve the drug-resistance profiles of HIV-1 non-nucleoside reverse transcriptase inhibitors (NNRTIs), a novel series of "dual-site" binding diarylpyrimidine (DAPY) derivatives targeting both the NNRTI adjacent site and NNRTIs binding pocket (NNIBP) were designed, synthesized, and evaluated for their anti-HIV potency in TZM-bl and MT-4 cells. Eight compounds exhibited moderate to excellent potencies in inhibiting wild-type (WT) HIV-1 replication with EC50 values ranging from 2.45 nM to 5.36 nM, and 14c (EC50 = 2.45 nM) proved to be the most promising inhibitor. Of note, 14c exhibited potent activity against the single mutant strain E138K (EC50 = 10.6 nM), being comparable with ETR (EC50 = 9.80 nM) and 3.5-fold more potent than that of compound 7 (EC50 = 37.3 nM). Moreover, 14c acted as a classical NNRTI with high affinity for WT HIV-1 RT (IC50 = 0.0589 µM). The detailed structure-activity relationships (SARs) of the representative compounds were also determined, and further supported by molecular dynamics simulation. Overall, we envision that the "dual-site"-binding NNRTIs have significant prospects and pave the way for the next round of rational design of potent anti-HIV-1 agents.

Comparison of Methods for Bulk Automated Simulation of Glycosidic Bond Conformations.

Six empirical force fields were tested for applicability to calculations for automated carbohydrate database filling. They were probed on eleven disaccharide molecules containing representative structural features from widespread classes of carbohydrates. The accuracy of each method was queried by predictions of nuclear Overhauser effects (NOEs) from conformational ensembles obtained from 50 to 100 ns molecular dynamics (MD) trajectories and their comparison to the published experimental data. Using various ranking schemes, it was concluded that explicit solvent MM3 MD yielded non-inferior NOE accuracy with newer GLYCAM-06, and ultimately PBE0-D3/def2-TZVP (Triple-Zeta Valence Polarized) Density Functional Theory (DFT) simulations. For seven of eleven molecules, at least one empirical force field with explicit solvent outperformed DFT in NOE prediction. The aggregate of characteristics (accuracy, speed, and compatibility) made MM3 dynamics with explicit solvent at 300 K the most favorable method for bulk generation of disaccharide conformation maps for massive database filling.

Estimating the Quality of 3D Protein Models Using the ModFOLD7 Server.

Assessing the accuracy of 3D models has become a keystone in the protein structure prediction field. ModFOLD7 is our leading resource for Estimates of Model Accuracy (EMA), which has been upgraded by integrating a number of the pioneering pure-single- and quasi-single-model approaches. Such an integration has given our latest version the strengths to accurately score and rank predicted models, with higher consistency compared to older EMA methods. Additionally, the server provides three options for producing global score estimates, depending on the requirements of the user: (1) ModFOLD7_rank, which is optimized for ranking/selection, (2) ModFOLD7_cor, which is optimized for correlations of predicted and observed scores, and (3) ModFOLD7 global for balanced performance. ModFOLD7 has been ranked among the top few EMA methods according to independent blind testing by the CASP13 assessors. Another evaluation resource for ModFOLD7 is the CAMEO project, where the method is continuously automatically evaluated, showing a significant improvement compared to our previous versions. The ModFOLD7 server is freely available at http://www.reading.ac.uk/bioinf/ModFOLD/ .

Molecular Dynamics Flexible Fitting: All You Want to Know About Resolution Exchange.

In recent years, owing to the advances in instrumentation, cryo-EM has emerged as the go-to tool for obtaining high-resolution structures of biomolecular systems. However, building three-dimensional atomic structures of biomolecules from these high-resolution maps remains a concern for the traditional map-guided structure-determination schemes. Recently, we developed a computational tool, Resolution Exchange Molecular Dynamics Flexible Fitting (ReMDFF) to address this problem by re-refining a search model against a series of maps of progressively higher resolutions, which ends with the original experimental resolution (Wang et al., J Struct Biol 204(2):319-328, 2018). In this chapter, we present a step-by-step outline for preparing, executing, and analyzing ReMDFF refinements of simple proteins and multimeric complexes. The structure determination of carbon monoxide dehydrogenase and Mg2+-channel CorA are employed as case studies. All scripts are provided via GitHub (Vant, Resolution exchange molecular dynamics flexible fitting (ReMDFF) all you want to know about flexible fitting, 2019, https://github.com/jvant/ReMDFF_Singharoy_Group.git ).

Protocols for Fast Simulations of Protein Structure Flexibility Using CABS-Flex and SURPASS.

Conformational flexibility of protein structures can play an important role in protein function. The flexibility is often studied using computational methods since experimental characterization can be difficult. Depending on protein system size, computational tools may require large computational resources or significant simplifications in the modeled systems to speed up calculations. In this work, we present the protocols for efficient simulations of flexibility of folded protein structures that use coarse-grained simulation tools of different resolutions: medium, represented by CABS-flex, and low, represented by SUPRASS. We test the protocols using a set of 140 globular proteins and compare the results with structure fluctuations observed in MD simulations, ENM modeling, and NMR ensembles. As demonstrated, CABS-flex predictions show high correlation to experimental and MD simulation data, while SURPASS is less accurate but promising in terms of future developments.

Protein Structure Prediction in CASP13 Using AWSEM-Suite.

Recently several techniques have emerged that significantly enhance the quality of predictions of protein tertiary structures. In this study, we describe the performance of AWSEM-Suite, an algorithm that incorporates template-based modeling and coevolutionary restraints with a realistic coarse-grained force field, AWSEM. With its roots in neural networks, AWSEM contains both physical and bioinformatical energies that have been optimized using energy landscape theory. AWSEM-Suite participated in CASP13 as a server predictor and generated reliable predictions for most targets. AWSEM-Suite ranked eighth in both the free-modeling category and the hard-to-model category and in one case provided the best submitted prediction. Here we critically discuss the prediction performance of AWSEM-Suite using several examples from different categories in CASP13. Structure prediction tests on these selected targets, two of them being hard-to-model targets, show that AWSEM-Suite can achieve high-resolution structure prediction after incorporating both template guidances and coevolutionary restraints even when homology is weak. For targets with reliable templates (template-easy category), introducing coevolutionary restraints sometimes damages the overall quality of the predictions. Free energy profile analyses demonstrate, however, that the incorporations of both of these evolutionarily informed terms effectively increase the funneling of the landscape toward native-like structures while still allowing sufficient flexibility to correct for discrepancies between the correct target structure and the provided guidance. In contrast to other predictors that are exclusively oriented toward structure prediction, the connection of AWSEM-Suite to a statistical mechanical basis and affiliated molecular dynamics and importance sampling simulations makes it suitable for functional explorations.

Analytical Method Using a Scaled Hypersphere Search for High-Dimensional Metadynamics Simulations.

Metadynamics (MTD) is one of the most effective methods for calculating the free energy surface and finding rare events. Nevertheless, numerous studies using MTD have been carried out using 3D or lower dimensional collective variables (CVs), as higher dimensional CVs require costly computational resources and the obtained results are too complex to understand the important events. The latter issue can be conveniently solved by utilizing the free energy reaction network (FERN), which is a graph structure consisting of edges of minimum free energy paths (MFEPs), nodes of equation (EQ) points, and transition state (TS) points. In the present article, a new method for exploring FERNs on high-dimensional CVs using MTD and the scaled hypersphere search (SHS) method is described. A test calculation based on the MTD-SHS simulation of met-enkephalin in explicit water with 7 CVs was conducted. As a result, 889 EQ points and 1805 TS points were found. The MTD-SHS approach can find MFEPs exhaustively; therefore, the FERNs can be estimated without any a priori knowledge of the EQ and TS points.

A New Lipid Force Field (FUJI).

To explore inhomogeneous and anisotropic systems such as lipid bilayers, the Lennard-Jones particle mesh Ewald (LJ-PME) method has been applied without a conventional isotropic dispersion correction. As the popular AMBER and CHARMM lipid force fields were developed using a cutoff scheme, their lipid bilayers unacceptably shrink when using the LJ-PME method. In this study, a new all-atom lipid force field (FUJI) was developed on the basis of the AMBER force-field scheme including the Lipid14 van der Waals parameters. Point charges were calculated using the restrained electrostatic potentials of many lipid conformers. Further, torsion energy profiles were calculated using high-level ab initio molecular orbitals (LCCSD(T)/Aug-cc-pVTZ//LMP2/Aug-cc-pVTZ), following which the molecular mechanical dihedral parameters were derived through a fast Fourier transform. By incorporation of these parameters into a new lipid force field without fitting experimental data, the desired lipid characteristics such as the area per lipid and lateral diffusion coefficients were obtained through GROMACS molecular dynamics simulations using the LJ-PME method and virtual hydrogen sites. The calculated area per lipid and lateral diffusion coefficients showed satisfactory agreement with experimental data. Furthermore, the electron-density profiles along the membrane normal were calculated for pure lipid bilayers, and the resulting membrane thicknesses agreed well with the experimental values. As the new lipid force field is compatible with FUJI for protein and small molecules, the new FUJI force field will offer accurate modeling for complex systems consisting of various membrane proteins and lipids.

Fine-Tuning of the AMBER RNA Force Field with a New Term Adjusting Interactions of Terminal Nucleotides.

Determination of RNA structural-dynamic properties is challenging for experimental methods. Thus, atomistic molecular dynamics (MD) simulations represent a helpful technique complementary to experiments. However, contemporary MD methods still suffer from limitations of force fields (ffs), including imbalances in the nonbonded ff terms. We have recently demonstrated that some improvement of state-of-the-art AMBER RNA ff can be achieved by adding a new term for H-bonding called gHBfix, which increases tuning flexibility and reduces risk of side-effects. Still, the first gHBfix version did not fully correct simulations of short RNA tetranucleotides (TNs). TNs are key benchmark systems due to availability of unique NMR data, although giving too much weight on improving TN simulations can easily lead to overfitting to A-form RNA. Here we combine the gHBfix version with another term called tHBfix, which separately treats H-bond interactions formed by terminal nucleotides. This allows to refine simulations of RNA TNs without affecting simulations of other RNAs. The approach is in line with adopted strategy of current RNA ffs, where the terminal nucleotides possess different parameters for terminal atoms than the internal nucleotides. Combination of gHBfix with tHBfix significantly improves the behavior of RNA TNs during well-converged enhanced-sampling simulations using replica exchange with solute tempering. TNs mostly populate canonical A-form like states while spurious intercalated structures are largely suppressed. Still, simulations of r(AAAA) and r(UUUU) TNs show some residual discrepancies with primary NMR data which suggests that future tuning of some other ff terms might be useful. Nevertheless, the tHBfix has a clear potential to improve modeling of key biochemical processes, where interactions of RNA single stranded ends are involved.

Semiclassical Vibrational Spectroscopy of Biological Molecules Using Force Fields.

Semiclassical spectroscopy is a practical way to get an accurately approximate quantum description of spectral features starting from ab initio molecular dynamics simulations. The computational bottleneck for the method is represented by the cost of ab initio potential, gradient, and Hessian matrix estimates. This drawback is particularly severe for biological systems due to their unique complexity and large dimensionality. The main goal of this manuscript is to demonstrate that quantum dynamics and spectroscopy, at the level of semiclassical approximation, are doable even for sizable biological systems. To this end, we investigate the possibility of performing semiclassical spectroscopy simulations when ab initio calculations are replaced by computationally cheaper force field evaluations. Both polarizable (AMOEBABIO18) and nonpolarizable (AMBER14SB) force fields are tested. Calculations of some particular vibrational frequencies of four nucleosides, i.e., uridine, thymidine, deoxyguanosine, and adenosine, show that ab initio simulations are accurate and widely applicable. Conversely, simulations based on AMBER14SB are limited to harmonic approximations, but those relying on AMOEBABIO18 yield acceptable semiclassical values if the investigated conformation has been included in the force field parametrization. The main conclusion is that AMOEBABIO18 may provide a viable route to assist semiclassical spectroscopy in the study of large biological molecules for which an ab initio approach is not computationally affordable.

Efficient Analytic Second Derivative of Electrostatic Embedding QM/MM Energy: Normal Mode Analysis of Plant Cryptochrome.

Analytic second derivatives of electrostatic embedding (EE) quantum mechanics/molecular mechanics (QM/MM) energy are important for performing vibrational analysis and simulating vibrational spectra of quantum systems interacting with an environment represented as a classical electrostatic potential. The main bottleneck of EE-QM/MM second derivatives is the solution of coupled perturbed equations for each MM atom perturbation. Here, we exploit the Q-vector method [J. Chem. Phys., 2019, 151, 041102] to workaround this bottleneck. We derive the full analytic second derivative of the EE-QM/MM energy, which allows us to compute QM, MM, and QM-MM Hessian blocks in an efficient and easy to implement manner. To show the capabilities of our method, we compute the normal modes for the full Arabidopsis thaliana plant cryptochrome. We show that the flavin adenine dinucleotide vibrations (QM subsystem) strongly mix with protein modes. We compute approximate vibronic couplings for the lowest bright transition, from which we extract spectral densities and the homogeneous broadening of FAD absorption spectrum in protein using vibrationally resolved electronic spectrum simulations.

Molecular Dynamics Simulations Reveal Canonical Conformations in Different pMHC/TCR Interactions.

The major defense system against microbial pathogens in vertebrates is the adaptive immune response and represents an effective mechanism in cancer surveillance. T cells represent an essential component of this complex system. They can recognize myriads of antigens as short peptides (p) originated from the intracellular degradation of foreign proteins presented by major histocompatibility complex (MHC) proteins. The clonotypic T-cell antigen receptor (TCR) is specialized in recognizing pMHC and triggering T cells immune response. It is still unclear how TCR engagement to pMHC is translated into the intracellular signal that initiates T-cell immune response. Some work has suggested the possibility that pMHC binding induces in the TCR conformational changes transmitted to its companion CD3 subunits that govern signaling. The conformational changes would promote phosphorylation of the CD3 complex ζ chain that initiates signal propagation intracellularly. Here, we used all-atom molecular dynamics simulations (MDs) of 500 ns to analyze the conformational behavior of three TCRs (1G4, ILA1 and ILA1α1ß1) interacting with the same MHC class I (HLA-A*02:01) bound to different peptides, and modelled in the presence of a lipid bilayer. Our data suggest a correlation between the conformations explored by the ß-chain constant regions and the T-cell response experimentally determined. In particular, independently by the TCR type involved in the interaction, the TCR activation seems to be linked to a specific zone of the conformational space explored by the ß-chain constant region. Moreover, TCR ligation restricts the conformational space the MHC class I groove.

Binding Dynamics of siRNA with Selected Lipopeptides: A Computer-Aided Study of the Effect of Lipopeptides' Functional Groups and Stereoisomerism.

The engineering issues pertaining to nanoparticle systems toward targeted gene therapies have not been fully probed. Recent experiments have identified specific structural characteristics of a novel class of lipopeptides (LP) that may lead to potent nanocarriers intended as RNAi therapeutics, albeit the molecular mechanism that underlies their performance remains unexplored. We conducted molecular dynamics simulations in atomistic detail coupled with free energy computations to study the dynamics and thermodynamics of an acrylate- and an epoxide-derived LP, members of the aforesaid class, upon their binding to siRNA in aqueous solution aiming at examining structure-potency relations. We found that the entropic part of the free energy of binding predominates; moreover, the first LP class tends to disrupt the Watson-Crick base pairing of siRNA, whereas the latter leaves the double helix intact. Moreover, the identified tug-of-war effect between LP-water and LP-siRNA hydrogen bonding in the supramolecular complex can underpin synthesis routes toward tuning the association dynamics. Our simulations on two diastereomers of the epoxide-derived LP showed significant structural and energetics differences upon binding, as a result of steric effects imposed by the different absolute configurations at their chiral centers. These findings may serve as crucial design parameters toward modulating the interplay between complex stability and ease of releasing the nucleic acid drug into the cell.

Markov State Model Analysis of Haloperidol Binding to the D3 Dopamine Receptor.

We have developed Markov state models (MSMs) and hidden Markov models (HMMs) that describe the binding of haloperidol to the D3 dopamine receptor. Haloperidol is an antipsychotic drug that binds with nanomolar affinity to the D3 dopamine receptor, where it functions as an inverse agonist. The models were constructed using an adaptive sampling approach from 519 individual molecular dynamics simulations totaling 122 µs of simulated time and encompass the entire drug binding process. They reveal short-lived metastable bound states and two distinct long-lived bound conformations that cannot be separated in affinity using our current methodology. This work extends the use of MSMs and HMMs to study ligand binding, which thus far has been limited to simpler systems.

Predicting Reactive Cysteines with Implicit-Solvent-Based Continuous Constant pH Molecular Dynamics in Amber.

Cysteines existing in the deprotonated thiolate form or having a tendency to become deprotonated are important players in enzymatic and cellular redox functions and frequently exploited in covalent drug design; however, most computational studies assume cysteines as protonated. Thus, developing an efficient tool that can make accurate and reliable predictions of cysteine protonation states is timely needed. We recently implemented a generalized Born (GB) based continuous constant pH molecular dynamics (CpHMD) method in Amber for protein pKa calculations on CPUs and GPUs. Here we benchmark the performance of GB-CpHMD for predictions of cysteine pKa's and reactivities using a data set of 24 proteins with both down- and upshifted cysteine pKa's. We found that 10 ns single-pH or 4 ns replica-exchange CpHMD titrations gave root-mean-square errors of 1.2-1.3 and correlation coefficients of 0.8-0.9 with respect to experiment. The accuracy of predicting thiolates or reactive cysteines at physiological pH with single-pH titrations is 86 or 81% with a precision of 100 or 90%, respectively. This performance well surpasses the traditional structure-based methods, particularly a widely used empirical pKa tool that gives an accuracy less than 50%. We discuss simulation convergence, dependence on starting structures, common determinants of the pKa downshifts and upshifts, and the origin of the discrepancies from the structure-based calculations. Our work suggests that CpHMD titrations can be performed on a desktop computer equipped with a single GPU card to predict cysteine protonation states for a variety of applications, from understanding biological functions to covalent drug design.