Time Scale Molecular Dynamics Research Articles

Structural insights into the activation of blood coagulation factor XI zymogen by thrombin: A computational molecular dynamics study

Activation of human blood coagulation factor XI zymogen to factor XIa plays a significant role in the upstream coagulation pathway, in which factor XIa activates factor IX zymogen. The mechanistic details of the proteolytic activation of factor XI by the activating enzyme thrombin are not well-understood at atomic level. In this study, we employed a combination of molecular docking and microsecond time-scale molecular dynamics simulations to identify the key regions of interaction between fXI and thrombin. The activating complex between the substrate and enzyme was modeled to represent the initial acylation step of the serine-protease hydrolysis mechanism. The proposed solution structural complex, fIX:fIIa, obtained from 3 microseconds of MD refinement, suggests that the activation of factor XI is mediated by thrombin's anion binding exosite-II interactions with A3 and A4 domains. We predict that the two positively charged arginine residues (Arg409 and Arg413) in the exosite-2 region, the β- and γ-insertion loops of thrombin play an important structural role in the initial activating complex between fXI and thrombin.

Structural Elucidation of Inter-CARD Interfaces involved in NOD2 Tandem CARD Association and RIP2 Recognition.

Nucleotide-binding and oligomerization domain-containing protein 2 (NOD2) recognizes the muramyl dipeptide and activates the NF-κB signaling cascade following its interaction with receptor-interacting protein 2 (RIP2) via caspase recruitment domains (CARDs). The NOD2-RIP2 interaction is not understood well due to inadequate structural information. Using comparative modeling and multimicrosecond timescale molecular dynamics simulations, we have demonstrated the association of NOD2-CARDs (CARDa-CARDb) and their interaction with RIP2CARD. Our results suggest that a negatively charged interface of NOD2CARDa and positively charged type-Ia interface of NOD2CARDb are crucial for CARDa-CARDb association and the type-Ia interface of NOD2CARDa and type-Ib interface of RIP2CARD predicted to be involved in 1:1 CARD-CARD interaction. Moreover, the direct interaction of NOD2CARDb with RIP2CARD signifies the importance of both CARDs of NOD2 in RIP2-mediated CARD-CARD interaction. Altogether, the structural results could help in understanding the underlying molecular details of the NOD2-RIP2 association in higher and lower eukaryotes.

Insight of displacement cascade evolution in gallium arsenide through molecular dynamics simulations

Displacement damage induced by space radiation is still a considerable challenge of expanding application of GaAs-based devices (e.g. solar cell, et al.) in space missions. In the current work, molecular dynamics (MD) simulations are applied to investigate displacement cascades evolution in GaAs at room temperature. A set of displacement events initiated by different energies (up to 50 keV) of primary knock-on atom (PKA) were simulated. An improved Wigner-Seitz (W-S) cell method was developed by introducing a simple criterion for more accurate identification of interstitials and antisites. Based on the present defect analysis method, the generation, evolution and distribution of point defects were discussed at first. Then, the mechanisms of both vacancy and interstitial cluster evolutions along simulation time were presented. The results indicate that both Ga and As type defects play the same crucial role among interstitials, antisites and vacancies. As the energy of PKA increases, the start points of heat spike phase postpone slightly. Moreover, the non-linear increase of surviving defects accompanied with the rise of the energy of PKA is considered as contribution of directly amorphous pockets. During the stage Ⅰ of cascade, as time progresses, different size clusters show up in succession from small size (i.e., between 2 and 5) to larger size. The medium and large vacancy clusters virtually only exist in the middle of cascades, and the majority of vacancies belong to the categorizes of single vacancy and small clusters, especially at the end of cascades. However, the medium and large interstitial clusters incline to hold up or grow continually at the MD time scale once they build up. Meanwhile, interstitials in larger clusters (30 ∼ ) account for more fraction as the energy of particles increases.

Multiple timescale molecular dynamics with very large time steps: avoidance of resonances

Multiple timescale molecular dynamics with very large time steps: avoidance of resonances

The Role of Grain Boundary Diffusion in the Solute Drag Effect.

Molecular dynamics (MD) simulations are applied to study solute drag by curvature-driven grain boundaries (GBs) in Cu–Ag solid solution. Although lattice diffusion is frozen on the MD timescale, the GB significantly accelerates the solute diffusion and alters the state of short-range order in lattice regions swept by its motion. The accelerated diffusion produces a nonuniform redistribution of the solute atoms in the form of GB clusters enhancing the solute drag by the Zener pinning mechanism. This finding points to an important role of lateral GB diffusion in the solute drag effect. A 1.5 at.%Ag alloying reduces the GB free energy by 10–20% while reducing the GB mobility coefficients by more than an order of magnitude. Given the greater impact of alloying on the GB mobility than on the capillary driving force, kinetic stabilization of nanomaterials against grain growth is likely to be more effective than thermodynamic stabilization aiming to reduce the GB free energy.

Thermal and Mechanochemical Tuning of the Porphyrin Singlet-Triplet Gap for Selective Energy Transfer Processes: A Molecular Dynamics Approach.

Molecular dynamics simulations provide fundamental knowledge on the reaction mechanism of a given simulated molecular process. Nevertheless, other methodologies based on the “static” exploration of potential energy surfaces are usually employed to firmly provide the reaction coordinate directly related to the reaction mechanism, as is the case in intrinsic reaction coordinates for thermally activated reactions. Photoinduced processes in molecular systems can also be studied with these two strategies, as is the case in the triplet energy transfer process. Triplet energy transfer is a fundamental photophysical process in photochemistry and photobiology, being for instance involved in photodynamic therapy, when generating the highly reactive singlet oxygen species. Here, we study the triplet energy transfer process between porphyrin, a prototypical energy transfer donor, and different biologically relevant acceptors, including molecular oxygen, carotenoids, and rhodopsin. The results obtained by means of nanosecond time-scale molecular dynamics simulations are compared to the “static” determination of the reaction coordinate for such a thermal process, leading to the distortions determining an effective energy transfer. This knowledge was finally applied to propose porphyrin derivatives for producing the required structural modifications in order to tune their singlet-triplet energy gap, thus introducing a mechanochemical description of the mechanism.

Structural origin of thermal shrinkage in soda-lime silicate glass below the glass transition temperature: A theoretical investigation by microsecond timescale molecular dynamics simulations.

Microscopic dynamical features in the relaxation of glass structures are one of the most important unsolved problems in condensed matter physics. Although the structural relaxation processes in the vicinity of glass transition temperature are phenomenologically expressed by the Kohlrausch-Williams-Watts function and the relaxation time can be successfully interpreted by Adam-Gibbs theory and/or Narayanaswamy's model, the atomic rearrangement, which is the origin of the volume change, and its driving force have not been elucidated. Using the microsecond time-scale molecular dynamics simulations, this study provides insights to quantitatively determine the origin of the thermal shrinkage below Tg in a soda-lime silicate glass. We found that during annealing below Tg, Na ions penetrate into the six-membered silicate rings, which remedies the acute O-O-O angles of the energetically unstable rings. The ring structure change makes the space to possess the cation inside the rings, but the ring volume is eventually reduced, which results in thermal shrinkage of the soda-lime silica glass. In conclusion, the dynamical structural relaxation due to the cation displacement evokes the overall volume relaxation at low temperature in the glassy material.

Addressing a Trapped High-Energy Water: Design and Synthesis of Highly Potent Pyrimidoindole-Based Glycogen Synthase Kinase-3β Inhibitors.

In small molecule binding, water is not a passive bystander but rather takes an active role in the binding site, which may be decisive for the potency of the inhibitor. Here, by addressing a high-energy water, we improved the IC50 value of our co-crystallized glycogen synthase kinase-3β (GSK-3β) inhibitor by nearly two orders of magnitude. Surprisingly, our results demonstrate that this high-energy water was not displaced by our potent inhibitor (S)-3-(3-((7-ethynyl-9H-pyrimido[4,5-b]indol-4-yl)(methyl)amino)piperidin-1-yl)propanenitrile ((S)-15, IC50 value of 6 nM). Instead, only a subtle shift in the location of this water molecule resulted in a dramatic decrease in the energy of this high-energy hydration site, as shown by the WaterMap analysis combined with microsecond timescale molecular dynamics simulations. (S)-15 demonstrated both a favorable kinome selectivity profile and target engagement in a cellular environment and reduced GSK-3 autophosphorylation in neuronal SH-SY5Y cells. Overall, our findings highlight that even a slight adjustment in the location of a high-energy water can be decisive for ligand binding.

Atomistic basis of force generation, translocation, and coordination in a viral genome packaging motor.

Double-stranded DNA viruses package their genomes into pre-assembled capsids using virally-encoded ASCE ATPase ring motors. We present the first atomic-resolution crystal structure of a multimeric ring form of a viral dsDNA packaging motor, the ATPase of the asccφ28 phage, and characterize its atomic-level dynamics via long timescale molecular dynamics simulations. Based on these results, and previous single-molecule data and cryo-EM reconstruction of the homologous φ29 motor, we propose an overall packaging model that is driven by helical-to-planar transitions of the ring motor. These transitions are coordinated by inter-subunit interactions that regulate catalytic and force-generating events. Stepwise ATP binding to individual subunits increase their affinity for the helical DNA phosphate backbone, resulting in distortion away from the planar ring towards a helical configuration, inducing mechanical strain. Subsequent sequential hydrolysis events alleviate the accumulated mechanical strain, allowing a stepwise return of the motor to the planar conformation, translocating DNA in the process. This type of helical-to-planar mechanism could serve as a general framework for ring ATPases.

Accelerated molecular dynamics simulation of vacancy diffusion in substitutional alloy with collective variable-driven hyperdynamics

It is not straightforward to investigate diffusion-limited phenomena in solid system by molecular dynamics (MD) simulation due to the limitation in time scale. Although many methods have been developed to extend the MD time scale by accelerating infrequent events, there is no established technique to treat diffusion in substitutional solid solution up to now since motions of solute atom, matrix atoms and vacancy interact each other complexly. In this study, an acceleration technique for the MD simulation of diffusion in substitutional solid solution is developed in conjunction with collective variable-driven hyperdynamics (CVHD). This novel technique is applied for the diffusion in Fe-Cr alloy as a representative substitutional alloy. Activation energy of vacancy diffusion from accelerated CVHD simulation basically agrees with that derived from MD simulation without acceleration, which confirms that the accelerated CVHD simulation reproduces the dynamics of vacancy diffusion in substitutional solid solution correctly.

Hamiltonian based resonance-free approach for enabling very large time steps in multiple time-scale molecular dynamics

Extended phase-space isokinetic methods in their deterministic [Minary et al., Phys. Rev. Lett. 93, 150201 (2004)] and stochastic forms [Leimkuhler et al., Mol. Phys. 111, 3579 (2013)] have proved tremendously successful in allowing multiple time-scale molecular dynamics simulations to be performed with very large time steps. These methods work by coupling the physical degrees of freedom to a set of Nosé-Hoover chain or Nosé-Hoover Langevin thermostats via an isokinetic constraint, which has the effect of avoiding resonance artifacts that plague multiple time-step algorithms. In this paper, we introduce a new resonance-free approach that achieves the same gains in time step but without the imposition of isokinetic constraints or the introduction of extended phase-space variables. Rather, we modify the physical Hamiltonian that effects the same regulation of resonances achieved by the isokinetic constraints. In so doing, we show that sampling errors can be controlled and performance improvements are possible within a simpler Hamiltonian framework. The method is demonstrated in simulations of the structure of liquid water and, in conjunction with enhanced sampling, in generation of the Ramachandran free-energy surface of the solvated alanine dipeptide.

Open Catalyst 2020 (OC20) Dataset and Community Challenges

Catalyst discovery and optimization is key to solving many societal and energy challenges including solar fuel synthesis, long-term energy storage, and renewable fertilizer production. Despite considerable effort by the catalysis community to apply machine learning models to the computational catalyst discovery process, it remains an open challenge to build models that can generalize across both elemental compositions of surfaces and adsorbate identity/configurations, perhaps because datasets have been smaller in catalysis than in related fields. To address this, we developed the OC20 dataset, consisting of 1,281,040 density functional theory (DFT) relaxations (∼264,890,000 single-point evaluations) across a wide swath of materials, surfaces, and adsorbates (nitrogen, carbon, and oxygen chemistries). We supplemented this dataset with randomly perturbed structures, short timescale molecular dynamics, and electronic structure analyses. The dataset comprises three central tasks indicative of day-to-day catalyst modeling and comes with predefined train/validation/test splits to facilitate direct comparisons with future model development efforts. We applied three state-of-the-art graph neural network models (CGCNN, SchNet, and DimeNet++) to each of these tasks as baseline demonstrations for the community to build on. In almost every task, no upper limit on model size was identified, suggesting that even larger models are likely to improve on initial results. The dataset and baseline models are both provided as open resources as well as a public leader board to encourage community contributions to solve these important tasks.

Characterizing the function of domain linkers in regulating the dynamics of multi-domain fusion proteins by microsecond molecular dynamics simulations and artificial intelligence.

Multi-domain proteins are not only formed through natural evolution but can also be generated by recombinant DNA technology. Because many fusion proteins can enhance the selectivity of cell targeting, these artificially produced molecules, called multi-specific biologics, are promising drug candidates, especially for immunotherapy. Moreover, the rational design of domain linkers in fusion proteins is becoming an essential step toward a quantitative understanding of the dynamics in these biopharmaceutics. We developed a computational framework to characterize the impacts of peptide linkers on the dynamics of multi-specific biologics. Specifically, we first constructed a benchmark containing six types of linkers that represent various lengths and degrees of flexibility and used them to connect two natural proteins as a test system. We then projected the microsecond dynamics of these proteins generated from Anton onto a coarse-grained conformational space. We further analyzed the similarity of dynamics among different proteins in this low-dimensional space by a neural-network-based classification model. Finally, we applied hierarchical clustering to place linkers into different subgroups based on the classification results. The clustering results suggest that the length of linkers, which is used to spatially separate different functional modules, plays the most important role in regulating the dynamics of this fusion protein. Given the same number of amino acids, linker flexibility functions as a regulator of protein dynamics. In summary, we illustrated that a new computational strategy can be used to study the dynamics of multi-domain fusion proteins by a combination of long timescale molecular dynamics simulation, coarse-grained feature extraction, and artificial intelligence.

Beryllium-driven structural evolution at the divertor surface

Erosion of the beryllium first wall material in tokamak reactors has been shown to result in transport and deposition on the tungsten divertor. Experimental studies of beryllium implantation in tungsten indicate that mixed W–Be intermetallic deposits can form, which have lower melting temperatures than tungsten and can trap tritium at higher rates. To better understand the formation and growth rate of these intermetallics, cumulative molecular dynamics (MD) simulations of both high and low energy beryllium deposition in tungsten were performed. In both cases, a W–Be mixed material layer (MML) emerged at the surface within several nanoseconds, either through energetic implantation or a thermally-activated exchange mechanism, respectively. While some ordering of the material into intermetallics occurred, fully ordered structures did not emerge from the deposition simulations. Targeted MD simulations of the MML to further study the rate of Be diffusion and intermetallic growth rates indicate that for both cases, the gradual re-structuring of the material into an ordered intermetallic layer is beyond accessible MD time scales(⩽1 μs). However, the rapid formation of the MML within nanoseconds indicates that beryllium deposition can influence other plasma species interactions at the surface and begin to alter the tungsten material properties. Therefore, beryllium deposition on the divertor surface, even in small amounts, is likely to cause significant changes in plasma-surface interactions and will need to be considered in future studies.

Structural Changes in Kir3.2 Channels Leading to Loss of K+ Selectivity

Inwardly rectifying potassium (Kir) channels are essential for numerous physiological processes, including neuronal and cardiac excitability. Recently, the very rare disease Keppen-Lubinsky syndrome (KPLBS), caused by de novo heterozygous mutations in the Kir3.2 channel has been described. KPLBS leads to severe developmental and intellectual problems, microcephaly, tightly adherent skin and severe generalized lipodystrophy. Recent breakthroughs in x-ray crystallography and cryo-electron microscopy provide insights into the structure of inward rectifier channels, and enable detailed characterization of the structural effects of disease-causing mutations. Multi-µs time scale molecular dynamics simulations on the WT Kir3.2 and disease mutant channel provide insights into the structural changes, caused by the point mutation, and the molecular mechanisms, leading to loss of K+ selectivity.

Computational Study of PmHMGR Thiohemiacetal Breakdown Transition State

3-hydroxyl-3-methylglutaryl CoA reductase (HMGR) catalyzes the rate determining step in the isoprenoid pathway using two NADH molecules. Even though this enzyme is targeted by statins, the atomistic details of the mechanism are still unknown. This objective can be met by combination of computational biophysics and experimental validations through time-resolved Laue crystallography for all steps of the reaction. The catalysis of HMGR can be broken down into three major chemical transformations including two hydride transfers and a thiohemiacetal breakdown step. Moreover, since two cofactor molecules are required, a cofactor exchange process also occurs after the first hydride transfer and before the second hydride transfer. Thiohemiacetal is a known intermediate resulting from the first hydride transfer step. However, the product of the thiohemiacetal breakdown step hypothesized to be mevaldehyde, has not been experimentally validated. This theoretical study explores the atomic details of the enzyme and ligands for the thiohemiacetal breakdown step using the Pseudomonas mevalonii HMGR. Theozymes and QM/MM ONIOM models were used to identify the transition state were performed to reveal the catalytic residues. These models have revealed that the Nδ of His381 is involved in stabilizing the negative charge in the thiolate anion. Ser85 forms a hydrogen bonding network with the first amide of thiohemiacetal ligand, a function performed by His381 in the preceding step first hydride transfer. Transition state force fields (TSFF) using the quantum guided molecular mechanics (Q2MM) program will then be used to perform µsec timescale molecular dynamics simulations. This project focuses on not only developing computational biophysics methods for studying enzyme mechanisms but also uses time-resolved x-ray crystallography to obtain snapshots along the HMGR reaction pathway which involves multiple transformations of the ligand in the same active site, a large flap domain movement, and a cofactor exchange.

Omega-3 versus Omega-6 fatty acid availability is controlled by hydrophobic site geometries of phospholipase A2s

Human phospholipase A2s (PLA2) constitute a superfamily of enzymes that hydrolyze the sn-2 acyl-chain of glycerophospholipids, producing lysophospholipids and free fatty acids. Each PLA2 enzyme type contributes to specific biological functions based on its expression, subcellular localization, and substrate specificity. Among the PLA2 superfamily, the cytosolic cPLA2 enzymes, calcium-independent iPLA2 enzymes, and secreted sPLA2 enzymes are implicated in many diseases, but a central issue is the preference for double-bond positions in polyunsaturated fatty acids (PUFAs) occupying the sn-2 position of membrane phospholipids. We demonstrate that each PLA2 has a unique preference between the specific omega-3 fatty acids eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA) and the omega-6 arachidonic acid (AA), which are the precursors of most proinflammatory and anti-inflammatory or resolving eicosanoids and related oxylipins. Surprisingly, we discovered that human cPLA2 selectively prefers AA, whereas iPLA2 prefers EPA, and sPLA2 prefers DHA as substrate. We determined the optimal binding of each phospholipid substrate in the active site of each PLA2 to explain these specificities. To investigate this, we utilized recently developed lipidomics-based LC-MS/MS and GC/MS assays to determine the sn-2 acyl chain specificity in mixtures of phospholipids. We performed μs timescale molecular dynamics (MD) simulations to reveal unique active site properties, especially how the precise hydrophobic cavity accommodation of the sn-2 acyl chain contributes to the stability of substrate binding and the specificity of each PLA2 for AA, EPA, or DHA. This study provides the first comprehensive picture of the unique substrate selectivity of each PLA2 for omega-3 and omega-6 fatty acids.

Conformational Landscapes of Halohydrin Dehalogenases and Their Accessible Active Site Tunnels

Halohydrin dehalogenases (HHDH) are industrially relevant biocatalysts exhibiting a promiscuous epoxide-ring opening reactivity in the presence of small nucleophiles, thus giving access to novel carbon–carbon, carbon–oxygen, carbon–nitrogen, and carbon–sulfur bonds. Recently, the repertoire of HHDH has been expanded, providing access to some novel HHDH subclasses exhibiting a broader epoxide substrate scope. In this work, we develop a computational approach based on the application of linear and non-linear dimensionality reduction techniques to long time-scale Molecular Dynamics (MD) simulations to study the HHDH conformational landscapes. We couple the analysis of the conformational landscapes to CAVER calculations to assess their impact on the active site tunnels and potential ability towards bulky epoxide ring opening reaction. Our study indicates that the analyzed HHDHs subclasses share a common breathing motion of the halide binding pocket, but present large deviations in the loops adjacent to the active site pocket and N-terminal regions. Such conformational differences affect the available tunnels for epoxide binding to the active site. The superior activity of the HHDH G subclass towards bulkier substrates is explained by the additional structural elements delimiting the active site region, its rich conformational heterogeneity, and the substantially wider and frequently observed active site tunnels. This study therefore provides key information for HHDH promiscuity and engineering.

Molecular Dynamics with Very Large Time Steps for the Calculation of Solvation Free Energies.

In multiple time scale molecular dynamics, the use of isokinetic constraints along with massive thermostatting has enabled the adoption of very large integration steps, well beyond the limits imposed by resonance artifacts in standard algorithms. In this work, we present two new contributions to this topic. First, we investigate the velocity distribution and the temperature-kinetic energy relationship associated with the isokinetic Nosé-Hoover family of methods, showing how they depend on the number of thermostats attached to each atomic degree of freedom. Second, we investigate the performance of these methods in the calculation of solvation free energies, the determination of which is often key for understanding the partition of a chemical species among distinct environments. We show how one can extract this property from canonical (constant-NVT) simulations and compare the result to experimental data obtained at a specific pressure. Finally, we demonstrate that large time steps can, in fact, be used to improve the efficiency of these calculations and that attaching multiple thermostats per degree of freedom is beneficial for effectively exploring the configurational space of a molecular system.

Ligand binding and global adaptation of the GlnPQ substrate binding domain 2 revealed by molecular dynamics simulations.

Substrate‐binding domains (SBD) are important structural elements of substrate transporters mediating the transport of essential molecules across the cell membrane. The SBD2 domain of the glutamine (GLN) transporter from bacteria consists of two domains D1 and D2 that bind GLN in the space between the domains in a closed conformation. In the absence of ligand, SBD2 adopts an open conformation with increased domain distance. In molecular dynamics (MD) simulations in the absence of ligands, no closing of the open conformation was observed on the MD time scale. Addition of GLN resulted in several reversible binding and unbinding events of GLN at the binding site on the D1 domain but did not induce domain closing indicating that binding and global domain closing do not occur simultaneously. The SBD2 structure remained in a closed state when starting from the GLN‐bound closed crystal structure and opened quickly to reach the open state upon removal of the GLN ligand. Free energy simulations to induce opening to closing indicated a barrier for closing in the absence and presence of ligand and a significant penalty for closing without GLN in the binding pocket. Simulations of a Leu480Ala mutation also indicate that an interaction of a C‐terminal D1‐tail471‐484 with a D2‐helix418‐427 (not contacting the substrate‐binding region) plays a decisive role for controlling the barrier of conformational switching in the SBD2 protein. The results allow us to derive a model of the molecular mechanism of substrate binding to SBD2 and associated conformational changes.