Independent MD Simulations Research Articles

Flap Dynamics in Pepsin-Like Aspartic Proteases: A Computational Perspective Using Plasmepsin-II and BACE-1 as Model Systems.

The flexibility of β hairpin structure known as the flap plays a key role in catalytic activity and substrate intake in pepsin-like aspartic proteases. Most of these enzymes share structural and sequential similarity. In this study, we have used apo Plm-II and BACE-1 as model systems. In the apo form of the proteases, a conserved tyrosine residue in the flap region remains in a dynamic equilibrium between the normal and flipped states through rotation of the χ1 and χ2 angles. Independent MD simulations of Plm-II and BACE-1 remained stuck either in the normal or flipped state. Metadynamics simulations using side-chain torsion angles (χ1 and χ2 of tyrosine) as collective variables sampled the transition between the normal and flipped states. Qualitatively, the two states were predicted to be equally populated. The normal and flipped states were stabilized by H-bond interactions to a tryptophan residue and to the catalytic aspartate, respectively. Further, mutation of tyrosine to an amino-acid with smaller side-chain, such as alanine, reduced the flexibility of the flap and resulted in a flap collapse (flap loses flexibility and remains stuck in a particular state). This is in accordance with previous experimental studies, which showed that mutation to alanine resulted in loss of activity in pepsin-like aspartic proteases. Our results suggest that the ring flipping associated with the tyrosine side-chain is the key order parameter that governs flap dynamics and opening of the binding pocket in most pepsin-like aspartic proteases.

Early stages of misfolding of PAP248-286 at two different pH values: An insight from molecular dynamics simulations

PAP248-286 peptides, which are highly abundant in human semen, aggregate and form amyloid fibrils that enhance HIV infection. Previous experimental studies have shown that the infection-promoting activity of PAP248-286 begins to increase well before amyloid formation takes place and that pH plays a key role in the enhancement of PAP248-286-related infection. Hence, understanding the early stages of misfolding of the PAP2482-86 peptide is crucial. To this end, we have performed 60 independent MD simulations for a total of 24 µs at two different pH values (4.2 and 7.2). Our data shows that early stages of misfolding of the PAP248-286 peptide is a multistage process and that the first step of the process is a transition from an “I-shaped” structure to a “U-shaped” structure. We further observed that the structure of PAP248-286 at the two different pH values shows significantly different features. At pH 4.2, the peptide has less intra-molecular H-bonds and a reduced α-helical content than at pH 7.2. Moreover, differences in intra-peptide residues contacts are also observed at the two pH values. Finally, free energy landscape analysis shows that there are more local minima in the energy surface of the peptide at pH 7.2 than at pH 4.2. Overall, the present study elucidates the early stages of misfolding of the PAP248-286 peptide at the atomic level, thus possibly opening new avenues in structure-based drug discovery against HIV infection.

Importance of protein intrinsic conformational dynamics and transient nature of non-covalent interactions in ligand binding affinity

We have recently identified BEN1 as a protein interactor of seryl-tRNA synthetase (SerRS) from model plant Arabidopsis thaliana. BEN1 contains an NADP+ binding domain and possesses acidic N-terminal extension essential for interaction with A. thaliana SerRS. This extension, specific for BEN1 homologues from Brassicaceae family, is solvent-exposed and distant to the nucleotide-binding site. We prepared a truncated BEN1 variant ΔN17BEN1 lacking the first 17 amino acid of this N-terminal extension as well as full-length BEN1 to investigate how the truncation affects the binding affinity towards coenzyme NADP+. By performing microscale thermophoresis (MST) experiments we have shown that both BEN1 variants bind the NADP+ cofactor, however, truncated BEN1 showed 34-fold higher affinity towards NADP+ indicating that its core protein structure is not just preserved but it binds NADP+ even stronger. To further corroborate the obtained results, we opted for a computational approach based on classical molecular dynamics simulations of both complexes. Our results have shown that both truncated and intact BEN1 variants form the same number of interactions with the NADP+ cofactor; however, it was the interaction occupancy that was affected. Namely, three independent MD simulations showed that the ΔN17BEN1 variant in complex with NADP+ has significantly higher interaction occupancy thus binds NADP+ with more than one order of magnitude higher affinity. Contrary to our expectations, the truncation of this distant region that does not communicate with the nucleotide-binding site didn't result in the gain of interaction but affected the intrinsic conformational dynamics which in turn fine-tuned the binding affinity by increasing the interaction occupancy and strength of the key conserved cation-π interaction between Arg69 and adenine of NADP+ and hydrogen bond between Ser244 and phosphate of NADP+.

Theoretical Studies of the Self Cleavage Pistol Ribozyme Mechanism

Ribozymes are huge complex biological catalysts composed of a combination of RNA and proteins. Nevertheless, there is a reduced number of small ribozymes, the self-cleavage ribozymes, that are formed just by RNA and, apparently, they existed in cells of primitive biological systems. Unveiling the details of these “fossils” enzymes can contribute not only to the understanding of the origins of life but also to the development of new simplified artificial enzymes. A computational study of the reactivity of the pistol ribozyme carried out by means of classical MD simulations and QM/MM hybrid calculations is herein presented to clarify its catalytic mechanism. Analysis of the geometries along independent MD simulations with different protonation states of the active site basic species reveals that only the canonical system, with no additional protonation changes, renders reactive conformations. A change in the coordination sphere of the Mg2+ ion has been observed during the simulations, which allows proposing a mechanism to explain the unique mode of action of the pistol ribozyme by comparison with other ribozymes. The present results are at the center of the debate originated from recent experimental and theoretical studies on pistol ribozyme.

Estimates of Electrical Conductivity from Molecular Dynamics Simulations: How to Invest the Computational Effort.

Although the electrical conductivity of an electrolyte can be estimated from the molecular dynamics trajectory, it is often a challenging task because of the need to obtain a substantial amount of data to ensure sufficient averaging. Here, we present an analysis on the convergence of results with the number of simulated trajectories. A series of molecular dynamics simulations have been performed for a model electrolyte (NaCl in water) and the Einstein relation has been used to calculate the electrical conductivity. The standard deviation of the conductivity estimates is relatively large compared to the mean value, and it has been shown that the off-diagonal contributions to the collective displacement of ions are responsible for large deviations between systems. It has been found that about 40 independent MD simulations may be required to reduce the errors. A procedure to improve the final estimate of the conductivity has been proposed.

Prediction of ligand binding mode among multiple cross-docking poses by molecular dynamics simulations.

We propose a method to identify the correct binding mode of a ligand with a protein among multiple predicted docking poses. Our method consists of two steps. First, five independent MD simulations with different initial velocities are performed for each docking pose, in order to evaluate its stability. If the root-mean-square deviations (RMSDs) of heavy atoms from the docking pose are larger than a given threshold (2.0Å) in all five parallel runs, the pose is filtered out and discarded. Then, we perform accurate all-atom binding free energy calculations for the residual poses only. The pose with the lowest binding free energy is identified as the correct pose. As a test case, we applied our method to a previously built cross-docking test set, which included 104 complex systems. We found that the present method could successfully identify the correct ligand binding mode for 72% (75/104) of the complexes for current test set. The possible reasons for the failure of the method in the other cases were investigated in detail, to enable future improvements.

Edge expansion parallel cascade selection molecular dynamics simulation for investigating large-amplitude collective motions of proteins.

We propose edge expansion parallel cascade selection molecular dynamics (eePaCS-MD) as an efficient adaptive conformational sampling method to investigate the large-amplitude motions of proteins without prior knowledge of the conformational transitions. In this method, multiple independent MD simulations are iteratively conducted from initial structures randomly selected from the vertices of a multi-dimensional principal component subspace. This subspace is defined by an ensemble of protein conformations sampled during previous cycles of eePaCS-MD. The edges and vertices of the conformational subspace are determined by solving the "convex hull problem." The sampling efficiency of eePaCS-MD is achieved by intensively repeating MD simulations from the vertex structures, which increases the probability of rare event occurrence to explore new large-amplitude collective motions. The conformational sampling efficiency of eePaCS-MD was assessed by investigating the open-close transitions of glutamine binding protein, maltose/maltodextrin binding protein, and adenylate kinase and comparing the results to those obtained using related methods. In all cases, the open-close transitions were simulated in ∼10 ns of simulation time or less, offering 1-3 orders of magnitude shorter simulation time compared to conventional MD. Furthermore, we show that the combination of eePaCS-MD and accelerated MD can further enhance conformational sampling efficiency, which reduced the total computational cost of observing the open-close transitions by at most 36%.

Wetting of pristine and functionalized nanocrystalline cellulose

Molecular dynamic simulations were used to investigate the wetting behavior of nanocrystalline cellulose allomorphs. Three 40ns independent MD simulations using CHARMM force field gave droplet profiles, whose contact angles were obtained with LBADSA plugin for ImageJ. Surface hydrophilicity was related to hydroxyl group availability for hydrogen bonding. Oxidation of C6 hydroxyl groups significantly boosts surface hydrophilicity, and even low levels of modification (~7%) leads to contact angles close to zero degrees. These data can be useful for biotechnological applications.

Dissociation Process of a MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and the Markov State Model.

Recently, we efficiently generated dissociation pathways of a protein-ligand complex without applying force bias with parallel cascade selection molecular dynamics (PaCS-MD) and showed that PaCS-MD in combination with the Markov state model (MSM) yielded a binding free energy comparable to experimental values. In this work, we applied the same procedure to a complex of MDM2 protein and the transactivation domain of p53 protein (TAD-p53), the latter of which is known to be very flexible in the unbound state. Using 30 independent MD simulations in PaCS-MD, we successfully generated 25 dissociation pathways of the complex, which showed complete or partial unfolding of the helical region of TAD-p53 during the dissociation process within an average simulation time of 154.8 ± 46.4 ns. The standard binding free energy obtained in combination with one-dimensional-, three-dimensional (3D)- or Cα-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the center of mass position of TAD-p53 relative to MDM2, the dissociation rate constant was calculated, which was comparable to those measured experimentally. Cα-MSM based on all Cα coordinates of TAD-p53 reproduced the experimentally measured standard binding free energy, and dissociation and association rate constants. We conclude that the combination of PaCS-MD and MSM offers an efficient computational procedure to calculate binding free energies and kinetic rates.

Conformational states of Zika virus non-structural protein 3 determined by molecular dynamics simulations with small-angle X-Ray scattering data

Zika virus (ZIKV) has become a great public health emergency. Its non-structural protein 3 (NS3) is a key enzyme in viral replication and has been considered as a potential therapeutic target. A conformational characterization of ZIKV NS3 is critical for a comprehensive understanding of its molecular interactions and functions. However, the high conformational flexibility of solution NS3 obstacles the structural characterization of NS3 solely from the experimental observable that averages over its heterogeneous conformations. Here, we employed replica exchange with solute tempering (REST) method to simulate the di-domain protein ZIKV NS3. Three independent MD simulations identified a conserved conformational ensemble of NS3, consisting of a major conformational state and several minor states from compact to loose conformations. The major state agrees well with the scattering profile from small-angle X-ray scattering (SAXS) experiments. Moreover, the simulated ensemble is supported by a direct data-fitting result that requires both short- and long-range structural contacts to recover the experimental data. We discussed the interplay between simulation and experiment in ensemble construction of flexible biomolecules and shed light on the physically derived conformational ensembles.

Computational simulations assessment of mutations impact on streptokinase (SK) from a group G streptococci with enhanced activity – insights into the functional roles of structural dynamics flexibility of SK and stabilization of SK–μplasmin catalytic complex

Streptokinase (SK), a plasminogen activator (PA) that converts inactive plasminogen (Pg) to plasmin (Pm), is a protein secreted by groups A, C, and G streptococci (GAS, GCS, and GGS, respectively), with high sequence divergence and functional heterogeneity. While roles of some residual changes in altered SK functionality are shown, the underlying structural mechanisms are less known. Herein, using computational approaches, we analyzed the conformational basis for the increased activity of SK from a GGS (SKG132) isolate with four natural residual substitutions (Ile33Phe, Arg45Gln, Asn228Lys, Phe287Ile) compared to the standard GCS (SKC). Using the crystal structure of SK.Pm catalytic complex as main template SKC.μPm catalytic complex was modeled through homology modeling process and validated by several online validation servers. Subsequently, SKG132.μPm structure was constructed by altering the corresponding residual substitutions. Results of three independent MD simulations showed increased RMSF values for SKG132.μPm, indicating the enhanced structural flexibility compared to SKC.μPm, specially in 170 and 250 loops and three regions: R1 (149–161), R2 (182–215) and R3 (224–229). In parallel, the average number of Hydrogen bonds in 170 loop, R2 and R3 (especially for Asn228Lys) of SKG132 compared to that of the SKC was decreased. Accordingly, residue interaction networks (RINs) analyses indicated that Asn228Lys might induce more level of structural flexibility by generation of free Lys256, while Phe287Ile and Ile33Phe enhanced the stabilization of the SKG132.μPm catalytic complex. These results denoted the potential role of the optimal dynamic state and stabilized catalytic complex for increased PA potencies of SK as a thrombolytic drug.

Uncertainty Quantification in Alchemical Free Energy Methods.

Alchemical free energy methods have gained much importance recently from several reports of improved ligand–protein binding affinity predictions based on their implementation using molecular dynamics simulations. A large number of variants of such methods implementing different accelerated sampling techniques and free energy estimators are available, each claimed to be better than the others in its own way. However, the key features of reproducibility and quantification of associated uncertainties in such methods have barely been discussed. Here, we apply a systematic protocol for uncertainty quantification to a number of popular alchemical free energy methods, covering both absolute and relative free energy predictions. We show that a reliable measure of error estimation is provided by ensemble simulation—an ensemble of independent MD simulations—which applies irrespective of the free energy method. The need to use ensemble methods is fundamental and holds regardless of the duration of time of the molecular dynamics simulations performed.

Nontargeted Parallel Cascade Selection Molecular Dynamics for Enhancing the Conformational Sampling of Proteins.

Nontargeted parallel cascade selection molecular dynamics (nt-PaCS-MD) is proposed as an efficient conformational sampling method to enhance the conformational transitions of proteins, which is an extension of the original targeted PaCS-MD (t-PaCS-MD). The original PaCS-MD comprises cycles of (i) selection of initial structures for multiple independent MD simulations toward a predetermined target and (ii) conformational sampling by the independent MDs. In nt-PaCS-MD, structures that significantly deviate from an average are regarded as candidates that have high potential to address other metastable states and are chosen as the initial structures in the selection. To select significantly deviated structures, we examine the root-mean-square deviation (RMSD) of snapshots generated from the average structure based on Gram-Schmidt orthogonalization. nt-PaCS-MD was applied to the folding of the mini-protein chignolin in implicit solvent and to the open-closed conformational transitions of T4 lysozyme (T4L) and glutamine binding protein (QBP) in explicit solvent. We show that nt-PaCS-MD can reach chignolin's native state and can also cause the open-closed transition of T4L and QBP on a nanosecond time scale, which are very efficient in terms of conformational sampling and comparable to that with t-PaCS-MD.

Molecular Simulations of Solved Co-crystallized X-Ray Structures Identify Action Mechanisms of PDEδ Inhibitors

PDEδ is a small protein that binds and controls the trafficking of RAS subfamily proteins. Its inhibition protects initiation of RAS signaling, and it is one of the common targets considered for oncological drug development. In this study, we used solved x-ray structures of inhibitor-bound PDEδ targets to investigate mechanisms of action of six independent all-atom MD simulations. An analysis of atomic simulations combined with the molecular mechanic-Poisson-Boltzmann solvent accessible surface area/generalized Born solvent accessible surface area calculations led to the identification of action mechanisms for a panel of novel PDEδ inhibitors. To the best of our knowledge, this study is one of the first in silico investigations on co-crystallized PDEδ protein. A detailed atomic-scale understanding of the molecular mechanism of PDEδ inhibition may assist in the design of novel PDEδ inhibitors. One of the most common side effects for diverse small molecules/kinase inhibitors is their off-target interactions with cardiac ion channels and human-ether-a-go-go channel specifically. Thus, all of the studied PDEδ inhibitors are also screened in silico at the central cavities of hERG1 potassium channels.

Theoretical insight into azobis-(benzo-18-crown-6) ether combined with the alkaline earth metal cations

It is desirable to research the photoswitchable supramolecular catalysts by theoretical methods. The azobis-(benzo-18-crown-6) ether (butterfly-AZO for short) combined with the alkaline earth metal cations (Mg2+, Ca2+, Sr2+ and Ba2+) investigated by using both density functional theory (DFT) and reactive molecular dynamics (reactive MD) simulations in this article. DFT calculations demonstrated that butterfly-AZO⋯Ba2+ is a suitable photocontrol catalytic candidate and its isomerization is mainly controlled by the NN rotational mechanism. Furthermore, 100 independent reactive MD simulations demonstrated that 100% trans isomers of butterfly-AZO⋯Ba2+ translated to the cis forms, in good agreement with the experimental data of ∼95%.

Molecular Dynamics Simulations of 441 Two-Residue Peptides in Aqueous Solution: Conformational Preferences and Neighboring Residue Effects with the Amber ff99SB-ildn-NMR Force Field.

Understanding the intrinsic conformational preferences of amino acids and the extent to which they are modulated by neighboring residues is a key issue for developing predictive models of protein folding and stability. Here we present the results of 441 independent explicit-solvent MD simulations of all possible two-residue peptides that contain the 20 standard amino acids with histidine modeled in both its neutral and protonated states. (3)J(HNHα) coupling constants and δ(Hα) chemical shifts calculated from the MD simulations correlate quite well with recently published experimental measurements for a corresponding set of two-residue peptides. Neighboring residue effects (NREs) on the average (3)J(HNHα) and δ(Hα) values of adjacent residues are also reasonably well reproduced, with the large NREs exerted experimentally by aromatic residues, in particular, being accurately captured. NREs on the secondary structure preferences of adjacent amino acids have been computed and compared with corresponding effects observed in a coil library and the average β-turn preferences of all amino acid types have been determined. Finally, the intrinsic conformational preferences of histidine, and its NREs on the conformational preferences of adjacent residues, are both shown to be strongly affected by the protonation state of the imidazole ring.

Theoretical design of visible light driven azobenzene-based photo-switching molecules

The preparation of switchable azobenzene derivatives driven by visible light is desirable for applications in biomolecular systems. o-R-substituted 4,4′-diacetamidoazobenzene derivatives (RH, CH3, OCH3 or OH) were investigated by using both density functional theory (DFT) and reactive molecular dynamics simulations. DFT calculations demonstrated that the nonplanar azo trans geometric structure, which caused by bulky groups tetra substituted in the ortho-position, is the key factor to enable the trans→cis transition with visible light. Furthermore, 100 independent reactive MD simulations demonstrated that 71% trans isomers of tetra o-OCH3-substituted 4,4′-diacetamidoazobenzene translated to cis, in good agreement with the experimental data.

Holo‐BtuF stabilizes the open conformation of the vitamin B12 ABC transporter BtuCD

While there is evidence that other ABC transporters can tell between empty and loaded substrate binding protein, reconstitution experiments suggest otherwise for the Escherichia coli vitamin B12 importer BtuCD-F. Here, we address the question of BtuCD-F substrate sensitivity in a combined protein-protein docking and molecular dynamics simulation approach. Starting from the BtuCD and holo-BtuF crystal structures, we model two holo-BtuCD-F docking complexes differing by a 180 degrees orientation of BtuF. One of these is similar to the apo-BtuCD-F crystal structure. Both docking complexes were embedded in a lipid/water environment to investigate their dynamics and BtuCD's conformational response to the presence and absence of BtuF, vitamin B12, and Mg-ATP in a series of 28 independent MD simulations. We find holo-BtuF stabilizing the open conformation of BtuCD, whereas the transporter begins to close again when BtuF or vitamin B12 is removed-suggesting BtuCD-F is capable of substrate sensitivity. We identified BtuC transmembrane helices 3 and 5, the L-loops and the adjacent helices comprised of BtuC residues 170-180 as hotspots of conformational change. We propose the latter to act as substrate sensors. BtuF-Trp44 appears to act as a lid on the vitamin B12 binding cleft in BtuF X-ray structures and protrudes into the BtuCD transport channel in one of our simulations, which might represent an initial step in vitamin B12 uptake. On an average, we observe subunit motions where the nucleotide binding domains approach each other while the transmembrane domains display an opening trend toward the periplasm.

Dynamics of Ca2+-saturated calmodulin D129N mutant studied by multiple molecular dynamics simulations.

Fifteen independent 1-nsec MD simulations of fully solvated Ca(2+) saturated calmodulin (CaM) mutant D129N were performed from different initial conditions to provide a sufficient statistical basis to gauge the significance of observed dynamical properties. In all MD simulations the four Ca(2+) ions remained in their binding sites, and retained a single water ligand as observed in the crystal structure. The coordination of Ca(2+) ions in EF-hands I, II, and III was sevenfold. In EF-hand IV, which was perturbed by the mutation of a highly conserved Asp129, an anomalous eightfold Ca(2+) coordination was observed. The Ca(2+) binding loop in EF-hand II was observed to dynamically sample conformations related to the Ca(2+)-free form. Repeated MD simulations implicate two well-defined conformations of Ca(2+) binding loop II, whereas similar effect was not observed for loops I, III, and IV. In 8 out of 15 MD simulations Ca(2+) binding loop II adopted an alternative conformation in which the Thr62 >C=O group was displaced from the Ca(2+) coordination by a water molecule, resulting in the Ca(2+) ion ligated by two water molecules. The alternative conformation of the Ca(2+) binding loop II appears related to the "closed" state involved in conformational exchange previously detected by NMR in the N-terminal domain fragment of CaM and the C-terminal domain fragment of the mutant E140Q. MD simulations suggest that conformations involved in microsecond exchange exist partially preformed on the nanosecond time scale.